Coverage for python/lsst/drp/tasks/gbdesAstrometricFit.py: 11%
1020 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-21 02:05 -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 astropy.table import vstack
33from sklearn.cluster import AgglomerativeClustering
34from smatch.matcher import Matcher
36import lsst.afw.geom as afwgeom
37import lsst.afw.table
38import lsst.geom
39import lsst.pex.config as pexConfig
40import lsst.pipe.base as pipeBase
41import lsst.sphgeom
42from lsst.meas.algorithms import (
43 LoadReferenceObjectsConfig,
44 ReferenceObjectLoader,
45 ReferenceSourceSelectorTask,
46)
47from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry
49from .build_camera import BuildCameraFromAstrometryTask
51__all__ = [
52 "calculate_apparent_motion",
53 "GbdesAstrometricFitConnections",
54 "GbdesAstrometricFitConfig",
55 "GbdesAstrometricFitTask",
56 "GbdesAstrometricMultibandFitConnections",
57 "GbdesAstrometricMultibandFitTask",
58 "GbdesGlobalAstrometricFitConnections",
59 "GbdesGlobalAstrometricFitConfig",
60 "GbdesGlobalAstrometricFitTask",
61 "GbdesGlobalAstrometricMultibandFitConnections",
62 "GbdesGlobalAstrometricMultibandFitTask",
63]
66def calculate_apparent_motion(catalog, refEpoch):
67 """Calculate shift from reference epoch to the apparent observed position
68 at another date.
70 This function calculates the shift due to proper motion combined with the
71 apparent motion due to parallax. This is not used in the
72 `GbdesAstrometricFitTask` or related child tasks, but is useful for
73 assessing results.
75 Parameters
76 ----------
77 catalog : `astropy.table.Table`
78 Table containing position, proper motion, parallax, and epoch for each
79 source, labeled by columns 'ra', 'dec', 'pmRA', 'pmDec', 'parallax',
80 and 'MJD'.
81 refEpoch : `astropy.time.Time`
82 Epoch of the reference position.
84 Returns
85 -------
86 apparentMotionRA : `np.ndarray` [`astropy.units.Quantity`]
87 RA shift in degrees.
88 apparentMotionDec : `np.ndarray` [`astropy.units.Quantity`]
89 Dec shift in degrees.
90 """
91 ra_rad = catalog["ra"].to(u.rad)
92 dec_rad = catalog["dec"].to(u.rad)
94 dt = (catalog["MJD"] - refEpoch).to(u.yr)
95 properMotionRA = catalog["pmRA"].to(u.deg / u.yr) * dt
96 properMotionDec = catalog["pmDec"].to(u.deg / u.yr) * dt
98 sun = astropy.coordinates.get_body("sun", time=catalog["MJD"])
99 frame = astropy.coordinates.GeocentricTrueEcliptic(equinox=catalog["MJD"])
100 sunLongitudes = sun.transform_to(frame).lon.radian
102 # These equations for parallax come from Equations 5.2 in Van de Kamp's
103 # book Stellar Paths. They differ from the parallax calculated in gbdes by
104 # ~0.01 mas, which is acceptable for QA and plotting purposes.
105 parallaxFactorRA = np.cos(wcsfit.EclipticInclination) * np.cos(ra_rad) * np.sin(sunLongitudes) - np.sin(
106 ra_rad
107 ) * np.cos(sunLongitudes)
108 parallaxFactorDec = (
109 np.sin(wcsfit.EclipticInclination) * np.cos(dec_rad)
110 - np.cos(wcsfit.EclipticInclination) * np.sin(ra_rad) * np.sin(dec_rad)
111 ) * np.sin(sunLongitudes) - np.cos(ra_rad) * np.sin(dec_rad) * np.cos(sunLongitudes)
112 parallaxDegrees = catalog["parallax"].to(u.degree)
113 parallaxCorrectionRA = parallaxDegrees * parallaxFactorRA
114 parallaxCorrectionDec = parallaxDegrees * parallaxFactorDec
116 apparentMotionRA = properMotionRA + parallaxCorrectionRA
117 apparentMotionDec = properMotionDec + parallaxCorrectionDec
119 return apparentMotionRA, apparentMotionDec
122def _make_ref_covariance_matrix(
123 refCat, inputUnit=u.radian, outputCoordUnit=u.marcsec, outputPMUnit=u.marcsec, version=1
124):
125 """Make a covariance matrix for the reference catalog including proper
126 motion and parallax.
128 The output is flattened to one dimension to match the format expected by
129 `gbdes`.
131 Parameters
132 ----------
133 refCat : `lsst.afw.table.SimpleCatalog`
134 Catalog including proper motion and parallax measurements.
135 inputUnit : `astropy.unit.core.Unit`
136 Units of the input catalog
137 outputCoordUnit : `astropy.unit.core.Unit`
138 Units required for the coordinates in the covariance matrix. `gbdes`
139 expects milliarcseconds.
140 outputPMUnit : `astropy.unit.core.Unit`
141 Units required for the proper motion/parallax in the covariance matrix.
142 `gbdes` expects milliarcseconds.
143 version : `int`
144 Version of the reference catalog. Version 2 includes covariance
145 measurements.
146 Returns
147 -------
148 cov : `list` [`float`]
149 Flattened output covariance matrix.
150 """
151 cov = np.zeros((len(refCat), 25))
152 if version == 1:
153 # Here is the standard ordering of components in the cov matrix,
154 # to match the PM enumeration in C++ code of gbdes package's Match.
155 # Each tuple gives: the array holding the 1d error,
156 # the string in Gaia column names for this
157 # the ordering in the Gaia catalog
158 # and the ordering of the tuples is the order we want in our cov matrix
159 raErr = (refCat["coord_raErr"]).to(outputCoordUnit).to_value()
160 decErr = (refCat["coord_decErr"]).to(outputCoordUnit).to_value()
161 raPMErr = (refCat["pm_raErr"]).to(outputPMUnit).to_value()
162 decPMErr = (refCat["pm_decErr"]).to(outputPMUnit).to_value()
163 parallaxErr = (refCat["parallaxErr"]).to(outputPMUnit).to_value()
164 stdOrder = (
165 (raErr, "ra", 0),
166 (decErr, "dec", 1),
167 (raPMErr, "pmra", 3),
168 (decPMErr, "pmdec", 4),
169 (parallaxErr, "parallax", 2),
170 )
172 k = 0
173 for i, pr1 in enumerate(stdOrder):
174 for j, pr2 in enumerate(stdOrder):
175 if pr1[2] < pr2[2]:
176 cov[:, k] = 0
177 elif pr1[2] > pr2[2]:
178 cov[:, k] = 0
179 else:
180 # diagnonal element
181 cov[:, k] = pr1[0] * pr2[0]
182 k = k + 1
184 elif version == 2:
185 positionParameters = ["coord_ra", "coord_dec", "pm_ra", "pm_dec", "parallax"]
186 units = [outputCoordUnit, outputCoordUnit, outputPMUnit, outputPMUnit, outputPMUnit]
187 k = 0
188 for i, pi in enumerate(positionParameters):
189 for j, pj in enumerate(positionParameters):
190 if i == j:
191 cov[:, k] = ((refCat[f"{pi}Err"].value) ** 2 * inputUnit**2).to(units[j] * units[j]).value
192 elif i > j:
193 cov[:, k] = (refCat[f"{pj}_{pi}_Cov"].value * inputUnit**2).to_value(units[i] * units[j])
194 else:
195 cov[:, k] = (refCat[f"{pi}_{pj}_Cov"].value * inputUnit**2).to_value(units[i] * units[j])
196 k += 1
197 return cov
200def _nCoeffsFromDegree(degree):
201 """Get the number of coefficients for a polynomial of a certain degree with
202 two variables.
204 This uses the general formula that the number of coefficients for a
205 polynomial of degree d with n variables is (n + d) choose d, where in this
206 case n is fixed to 2.
208 Parameters
209 ----------
210 degree : `int`
211 Degree of the polynomial in question.
213 Returns
214 -------
215 nCoeffs : `int`
216 Number of coefficients for the polynomial in question.
217 """
218 nCoeffs = int((degree + 2) * (degree + 1) / 2)
219 return nCoeffs
222def _degreeFromNCoeffs(nCoeffs):
223 """Get the degree for a polynomial with two variables and a certain number
224 of coefficients.
226 This is done by applying the quadratic formula to the
227 formula for calculating the number of coefficients of the polynomial.
229 Parameters
230 ----------
231 nCoeffs : `int`
232 Number of coefficients for the polynomial in question.
234 Returns
235 -------
236 degree : `int`
237 Degree of the polynomial in question.
238 """
239 degree = int(-1.5 + 0.5 * (1 + 8 * nCoeffs) ** 0.5)
240 return degree
243def _convert_to_ast_polymap_coefficients(coefficients):
244 """Convert vector of polynomial coefficients from the format used in
245 `gbdes` into AST format (see Poly2d::vectorIndex(i, j) in
246 gbdes/gbutil/src/Poly2d.cpp). This assumes two input and two output
247 coordinates.
249 Parameters
250 ----------
251 coefficients : `list`
252 Coefficients of the polynomials.
253 degree : `int`
254 Degree of the polynomial.
256 Returns
257 -------
258 astPoly : `astshim.PolyMap`
259 Coefficients in AST polynomial format.
260 """
261 polyArray = np.zeros((len(coefficients), 4))
262 N = len(coefficients) / 2
263 degree = _degreeFromNCoeffs(N)
265 for outVar in [1, 2]:
266 for i in range(degree + 1):
267 for j in range(degree + 1):
268 if (i + j) > degree:
269 continue
270 vectorIndex = int(((i + j) * (i + j + 1)) / 2 + j + N * (outVar - 1))
271 polyArray[vectorIndex, 0] = coefficients[vectorIndex]
272 polyArray[vectorIndex, 1] = outVar
273 polyArray[vectorIndex, 2] = i
274 polyArray[vectorIndex, 3] = j
276 astPoly = astshim.PolyMap(polyArray, 2, options="IterInverse=1,NIterInverse=10,TolInverse=1e-7")
277 return astPoly
280def _get_instruments(inputVisitSummaries):
281 """Make `wcsfit.Instrument` objects for all of the instruments and filters
282 used for the input visits. This also returns the indices to match the
283 visits to the instrument/filter used.
285 Parameters
286 ----------
287 inputVisitSummaries: `list` [`lsst.afw.table.ExposureCatalog`]
288 List of catalogs with per-detector summary information.
290 Returns
291 -------
292 instruments : `list` [`wcsfit.Instrument`]
293 List of instrument objects.
294 instrumentIndices : `list` [`int`]
295 Indices matching each visit to the instrument/filter used.
296 """
297 instruments = []
298 filters = []
299 instrumentIndices = []
300 for visitSummary in inputVisitSummaries:
301 filter = visitSummary[0]["physical_filter"]
302 instrumentName = visitSummary[0].getVisitInfo().instrumentLabel
303 if filter not in filters:
304 filters.append(filter)
305 filter_instrument = wcsfit.Instrument(instrumentName)
306 filter_instrument.band = filter
307 instruments.append(filter_instrument)
308 instrumentIndices.append(filters.index(filter))
309 return instruments, instrumentIndices
312class CholeskyError(pipeBase.AlgorithmError):
313 """Raised if the Cholesky decomposition in the model fit fails."""
315 def __init__(self) -> None:
316 super().__init__(
317 "Cholesky decomposition failed, likely because data is not sufficient to constrain the model."
318 )
320 @property
321 def metadata(self) -> dict:
322 """There is no metadata associated with this error."""
323 return {}
326class GbdesAstrometricFitConnections(
327 pipeBase.PipelineTaskConnections,
328 dimensions=("skymap", "tract", "instrument", "physical_filter"),
329 defaultTemplates={
330 "outputName": "gbdesAstrometricFit",
331 },
332):
333 """Middleware input/output connections for task data."""
335 inputCatalogRefs = pipeBase.connectionTypes.Input(
336 doc="Source table in parquet format, per visit.",
337 name="preSourceTable_visit",
338 storageClass="DataFrame",
339 dimensions=("instrument", "visit"),
340 deferLoad=True,
341 multiple=True,
342 )
343 inputVisitSummaries = pipeBase.connectionTypes.Input(
344 doc=(
345 "Per-visit consolidated exposure metadata built from calexps. "
346 "These catalogs use detector id for the id and must be sorted for "
347 "fast lookups of a detector."
348 ),
349 name="visitSummary",
350 storageClass="ExposureCatalog",
351 dimensions=("instrument", "visit"),
352 multiple=True,
353 )
354 referenceCatalog = pipeBase.connectionTypes.PrerequisiteInput(
355 doc="The astrometry reference catalog to match to loaded input catalog sources.",
356 name="gaia_dr3_20230707",
357 storageClass="SimpleCatalog",
358 dimensions=("skypix",),
359 deferLoad=True,
360 multiple=True,
361 )
362 colorCatalog = pipeBase.connectionTypes.Input(
363 doc="The catalog of magnitudes to match to input sources.",
364 name="fgcm_Cycle4_StandardStars",
365 storageClass="SimpleCatalog",
366 dimensions=("instrument",),
367 )
368 inputCameraModel = pipeBase.connectionTypes.PrerequisiteInput(
369 doc="Camera parameters to use for 'device' part of model",
370 name="gbdesAstrometricFit_cameraModel",
371 storageClass="ArrowNumpyDict",
372 dimensions=("instrument", "physical_filter"),
373 )
374 inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
375 doc="Input camera to use when constructing camera from astrometric model.",
376 name="camera",
377 storageClass="Camera",
378 dimensions=("instrument",),
379 isCalibration=True,
380 )
381 outputWcs = pipeBase.connectionTypes.Output(
382 doc=(
383 "Per-tract, per-visit world coordinate systems derived from the fitted model."
384 " These catalogs only contain entries for detectors with an output, and use"
385 " the detector id for the catalog id, sorted on id for fast lookups of a detector."
386 ),
387 name="{outputName}SkyWcsCatalog",
388 storageClass="ExposureCatalog",
389 dimensions=("instrument", "visit", "skymap", "tract"),
390 multiple=True,
391 )
392 outputCatalog = pipeBase.connectionTypes.Output(
393 doc=(
394 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
395 "plane coordinates and chisq values."
396 ),
397 name="{outputName}_fitStars",
398 storageClass="ArrowNumpyDict",
399 dimensions=("instrument", "skymap", "tract", "physical_filter"),
400 )
401 starCatalog = pipeBase.connectionTypes.Output(
402 doc=(
403 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
404 "fitProperMotion is True."
405 ),
406 name="{outputName}_starCatalog",
407 storageClass="ArrowNumpyDict",
408 dimensions=("instrument", "skymap", "tract", "physical_filter"),
409 )
410 modelParams = pipeBase.connectionTypes.Output(
411 doc="WCS parameters and covariance.",
412 name="{outputName}_modelParams",
413 storageClass="ArrowNumpyDict",
414 dimensions=("instrument", "skymap", "tract", "physical_filter"),
415 )
416 outputCameraModel = pipeBase.connectionTypes.Output(
417 doc="Camera parameters to use for 'device' part of model",
418 name="{outputName}_cameraModel",
419 storageClass="ArrowNumpyDict",
420 dimensions=("instrument", "skymap", "tract", "physical_filter"),
421 )
422 camera = pipeBase.connectionTypes.Output(
423 doc="Camera object constructed using the per-detector part of the astrometric model",
424 name="{outputName}Camera",
425 storageClass="Camera",
426 dimensions=("instrument", "skymap", "tract", "physical_filter"),
427 )
428 dcrCoefficients = pipeBase.connectionTypes.Output(
429 doc="Per-visit coefficients for DCR correction.",
430 name="{outputName}_dcrCoefficients",
431 storageClass="ArrowNumpyDict",
432 dimensions=("instrument", "skymap", "tract", "physical_filter"),
433 )
435 def getSpatialBoundsConnections(self):
436 return ("inputVisitSummaries",)
438 def __init__(self, *, config=None):
439 super().__init__(config=config)
441 if self.config.healpix is not None:
442 self.dimensions.remove("tract")
443 self.dimensions.remove("skymap")
444 healpixName = f"healpix{self.config.healpix}"
445 self.dimensions.add(healpixName)
446 self.outputWcs = dataclasses.replace(
447 self.outputWcs, dimensions=("instrument", "visit", healpixName)
448 )
449 self.outputCatalog = dataclasses.replace(
450 self.outputCatalog, dimensions=("instrument", "physical_filter", healpixName)
451 )
452 self.starCatalog = dataclasses.replace(
453 self.starCatalog, dimensions=("instrument", "physical_filter", healpixName)
454 )
455 self.modelParams = dataclasses.replace(
456 self.modelParams, dimensions=("instrument", "physical_filter", healpixName)
457 )
458 self.outputCameraModel = dataclasses.replace(
459 self.outputCameraModel, dimensions=("instrument", "physical_filter", healpixName)
460 )
461 self.camera = dataclasses.replace(
462 self.camera, dimensions=("instrument", "physical_filter", healpixName)
463 )
464 self.dcrCoefficients = dataclasses.replace(
465 self.dcrCoefficients, dimensions=("instrument", "physical_filter", healpixName)
466 )
468 if not self.config.useColor:
469 self.inputs.remove("colorCatalog")
470 self.outputs.remove("dcrCoefficients")
471 if not self.config.saveModelParams:
472 self.outputs.remove("modelParams")
473 if not self.config.useInputCameraModel:
474 self.prerequisiteInputs.remove("inputCameraModel")
475 if not self.config.saveCameraModel:
476 self.outputs.remove("outputCameraModel")
477 if not self.config.saveCameraObject:
478 self.prerequisiteInputs.remove("inputCamera")
479 self.outputs.remove("camera")
482class GbdesAstrometricFitConfig(
483 pipeBase.PipelineTaskConfig, pipelineConnections=GbdesAstrometricFitConnections
484):
485 """Configuration for GbdesAstrometricFitTask"""
487 sourceSelector = sourceSelectorRegistry.makeField(
488 doc="How to select sources for cross-matching.", default="science"
489 )
490 referenceSelector = pexConfig.ConfigurableField(
491 target=ReferenceSourceSelectorTask,
492 doc="How to down-select the loaded astrometry reference catalog.",
493 )
494 referenceFilter = pexConfig.Field(
495 dtype=str,
496 doc="Name of filter to load from reference catalog. This is a required argument, although the values"
497 "returned are not used.",
498 default="phot_g_mean",
499 )
500 setRefEpoch = pexConfig.Field(
501 dtype=float,
502 doc="Set the reference epoch to a fixed value in MJD (if None, median observation date is used)",
503 default=None,
504 optional=True,
505 )
506 applyRefCatProperMotion = pexConfig.Field(
507 dtype=bool,
508 doc="Apply proper motion to shift reference catalog to epoch of observations.",
509 default=True,
510 )
511 matchRadius = pexConfig.Field(
512 doc="Matching tolerance between associated objects (arcseconds).", dtype=float, default=1.0
513 )
514 minMatches = pexConfig.Field(
515 doc="Number of matches required to keep a source object.", dtype=int, default=2
516 )
517 allowSelfMatches = pexConfig.Field(
518 doc="Allow multiple sources from the same visit to be associated with the same object.",
519 dtype=bool,
520 default=False,
521 )
522 sourceFluxType = pexConfig.Field(
523 dtype=str,
524 doc="Source flux field to use in source selection and to get fluxes from the catalog.",
525 default="apFlux_12_0",
526 )
527 systematicError = pexConfig.Field(
528 dtype=float,
529 doc=(
530 "Systematic error padding added in quadrature for the science catalogs (marcsec). The default"
531 "value is equivalent to 0.02 pixels for HSC."
532 ),
533 default=0.0034,
534 )
535 referenceSystematicError = pexConfig.Field(
536 dtype=float,
537 doc="Systematic error padding added in quadrature for the reference catalog (marcsec).",
538 default=0.0,
539 )
540 useColor = pexConfig.Field(
541 dtype=bool,
542 doc="Use color information to correct for differential chromatic refraction.",
543 default=False,
544 )
545 color = pexConfig.ListField(
546 dtype=str,
547 doc="The bands to use for calculating color.",
548 default=["g", "i"],
549 listCheck=(lambda x: (len(x) == 2) and (len(set(x)) == len(x))),
550 )
551 referenceColor = pexConfig.Field(
552 dtype=float,
553 doc="The color for which DCR is defined as zero.",
554 default=0.61,
555 )
556 modelComponents = pexConfig.ListField(
557 dtype=str,
558 doc=(
559 "List of mappings to apply to transform from pixels to sky, in order of their application."
560 "Supported options are 'INSTRUMENT/DEVICE' and 'EXPOSURE'."
561 ),
562 default=["INSTRUMENT/DEVICE", "EXPOSURE"],
563 )
564 deviceModel = pexConfig.ListField(
565 dtype=str,
566 doc=(
567 "List of mappings to apply to transform from detector pixels to intermediate frame. Map names"
568 "should match the format 'BAND/DEVICE/<map name>'."
569 ),
570 default=["BAND/DEVICE/poly"],
571 )
572 exposureModel = pexConfig.ListField(
573 dtype=str,
574 doc=(
575 "List of mappings to apply to transform from intermediate frame to sky coordinates. Map names"
576 "should match the format 'EXPOSURE/<map name>'."
577 ),
578 default=["EXPOSURE/poly"],
579 )
580 devicePolyOrder = pexConfig.Field(dtype=int, doc="Order of device polynomial model.", default=4)
581 exposurePolyOrder = pexConfig.Field(dtype=int, doc="Order of exposure polynomial model.", default=6)
582 fitProperMotion = pexConfig.Field(dtype=bool, doc="Fit the proper motions of the objects.", default=False)
583 excludeNonPMObjects = pexConfig.Field(
584 dtype=bool, doc="Exclude reference objects without proper motion/parallax information.", default=True
585 )
586 fitReserveFraction = pexConfig.Field(
587 dtype=float, default=0.2, doc="Fraction of objects to reserve from fit for validation."
588 )
589 fitReserveRandomSeed = pexConfig.Field(
590 dtype=int,
591 doc="Set the random seed for selecting data points to reserve from the fit for validation.",
592 default=1234,
593 )
594 saveModelParams = pexConfig.Field(
595 dtype=bool,
596 doc=(
597 "Save the parameters and covariance of the WCS model. Default to "
598 "false because this can be very large."
599 ),
600 default=False,
601 )
602 useInputCameraModel = pexConfig.Field(
603 dtype=bool,
604 doc=(
605 "Use a preexisting model for the 'device' part of the model. When true, the device part of the"
606 " model will be held fixed in the fitting process."
607 ),
608 default=False,
609 )
610 saveCameraModel = pexConfig.Field(
611 dtype=bool,
612 doc="Save the 'device' part of the model to be used as input in future runs.",
613 default=False,
614 )
615 buildCamera = pexConfig.ConfigurableField(
616 target=BuildCameraFromAstrometryTask, doc="Subtask to build camera from astrometric model."
617 )
618 saveCameraObject = pexConfig.Field(
619 dtype=bool,
620 doc="Build and output an lsst.afw.cameraGeom.Camera object using the fit per-detector model.",
621 default=False,
622 )
623 clipThresh = pexConfig.Field(
624 dtype=float,
625 doc="Threshold for clipping outliers in the fit (in standard deviations)",
626 default=5.0,
627 )
628 clipFraction = pexConfig.Field(
629 dtype=float,
630 doc="Minimum fraction of clipped sources that triggers a new fit iteration.",
631 default=0.0,
632 )
633 healpix = pexConfig.Field(
634 dtype=int,
635 doc="Run using all visits overlapping a healpix pixel with this order instead of a tract. Order 3 "
636 "corresponds to pixels with angular size of 7.329 degrees.",
637 optional=True,
638 default=None,
639 )
640 minDetectorFraction = pexConfig.Field(
641 dtype=float,
642 doc=(
643 "Minimum fraction of detectors with valid WCSs per visit required to include the visit in the "
644 "fit."
645 ),
646 default=0.25,
647 )
649 def setDefaults(self):
650 # Use only stars because aperture fluxes of galaxies are biased and
651 # depend on seeing.
652 self.sourceSelector["science"].doUnresolved = True
653 self.sourceSelector["science"].unresolved.name = "sizeExtendedness"
655 # Use only isolated sources.
656 self.sourceSelector["science"].doIsolated = True
657 self.sourceSelector["science"].isolated.parentName = "parentSourceId"
658 self.sourceSelector["science"].isolated.nChildName = "deblend_nChild"
659 # Do not use either flux or centroid measurements with flags,
660 # chosen from the usual QA flags for stars.
661 self.sourceSelector["science"].doFlags = True
662 badFlags = [
663 "pixelFlags_edge",
664 "pixelFlags_saturated",
665 "pixelFlags_interpolatedCenter",
666 "pixelFlags_interpolated",
667 "pixelFlags_crCenter",
668 "pixelFlags_bad",
669 "hsmPsfMoments_flag",
670 f"{self.sourceFluxType}_flag",
671 ]
672 self.sourceSelector["science"].flags.bad = badFlags
674 # Use only primary sources.
675 self.sourceSelector["science"].doRequirePrimary = True
677 self.sourceSelector["science"].doSignalToNoise = True
678 self.sourceSelector["science"].signalToNoise.minimum = 8.0
679 self.sourceSelector["science"].signalToNoise.maximum = 1000.0
680 self.sourceSelector["science"].signalToNoise.fluxField = self.sourceFluxType + "_instFlux"
681 self.sourceSelector["science"].signalToNoise.errField = self.sourceFluxType + "_instFluxErr"
683 def validate(self):
684 super().validate()
686 # Check if all components of the device and exposure models are
687 # supported.
688 for component in self.deviceModel:
689 mapping = component.split("/")[-1]
690 if mapping not in ["poly", "identity"]:
691 raise pexConfig.FieldValidationError(
692 GbdesAstrometricFitConfig.deviceModel,
693 self,
694 f"deviceModel component {component} is not supported.",
695 )
697 for component in self.exposureModel:
698 mapping = component.split("/")[-1]
699 if mapping not in ["poly", "identity", "dcr"]:
700 raise pexConfig.FieldValidationError(
701 GbdesAstrometricFitConfig.exposureModel,
702 self,
703 f"exposureModel component {component} is not supported.",
704 )
706 if self.saveCameraModel and self.useInputCameraModel:
707 raise pexConfig.FieldValidationError(
708 GbdesAstrometricFitConfig.saveCameraModel,
709 self,
710 "saveCameraModel and useInputCameraModel cannot both be true.",
711 )
713 if self.saveCameraObject and (self.exposurePolyOrder != 1):
714 raise pexConfig.FieldValidationError(
715 GbdesAstrometricFitConfig.saveCameraObject,
716 self,
717 "If saveCameraObject is True, exposurePolyOrder must be set to 1.",
718 )
721class GbdesAstrometricFitTask(pipeBase.PipelineTask):
722 """Calibrate the WCS across multiple visits of the same field using the
723 GBDES package.
724 """
726 ConfigClass = GbdesAstrometricFitConfig
727 _DefaultName = "gbdesAstrometricFit"
729 def __init__(self, **kwargs):
730 super().__init__(**kwargs)
731 self.makeSubtask("sourceSelector")
732 self.makeSubtask("referenceSelector")
733 if self.config.saveCameraObject:
734 self.makeSubtask("buildCamera")
736 def runQuantum(self, butlerQC, inputRefs, outputRefs):
737 # We override runQuantum to set up the refObjLoaders
738 inputs = butlerQC.get(inputRefs)
740 instrumentName = butlerQC.quantum.dataId["instrument"]
742 # Ensure the inputs are in a consistent and deterministic order
743 inputCatVisits = np.array([inputCat.dataId["visit"] for inputCat in inputs["inputCatalogRefs"]])
744 inputs["inputCatalogRefs"] = [inputs["inputCatalogRefs"][v] for v in inputCatVisits.argsort()]
745 inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]])
746 inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()]
747 inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog])
748 inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()]
749 inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]])
750 inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()]
752 refConfig = LoadReferenceObjectsConfig()
753 if self.config.applyRefCatProperMotion:
754 refConfig.requireProperMotion = True
755 refObjectLoader = ReferenceObjectLoader(
756 dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs],
757 refCats=inputs.pop("referenceCatalog"),
758 config=refConfig,
759 log=self.log,
760 )
762 nCores = butlerQC.resources.num_cores
763 self.log.info("Running with nCores = %d", nCores)
765 if self.config.useColor:
766 colorCatalog = inputs.pop("colorCatalog")
767 else:
768 colorCatalog = None
769 try:
770 output = self.run(
771 **inputs,
772 instrumentName=instrumentName,
773 refObjectLoader=refObjectLoader,
774 colorCatalog=colorCatalog,
775 nCores=nCores,
776 )
777 except pipeBase.AlgorithmError as e:
778 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
779 # No partial outputs for butler to put
780 raise error from e
782 wcsOutputRefDict = {outWcsRef.dataId["visit"]: outWcsRef for outWcsRef in outputRefs.outputWcs}
783 for visit, outputWcs in output.outputWcss.items():
784 butlerQC.put(outputWcs, wcsOutputRefDict[visit])
785 butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
786 butlerQC.put(output.starCatalog, outputRefs.starCatalog)
787 if self.config.saveModelParams:
788 butlerQC.put(output.modelParams, outputRefs.modelParams)
789 if self.config.saveCameraModel:
790 butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel)
791 if self.config.saveCameraObject:
792 butlerQC.put(output.camera, outputRefs.camera)
793 if self.config.useColor:
794 butlerQC.put(output.colorParams, outputRefs.dcrCoefficients)
795 if output.partialOutputs:
796 e = RuntimeError("Some visits were dropped because data was insufficient to fit model.")
797 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
798 raise error from e
800 def run(
801 self,
802 inputCatalogRefs,
803 inputVisitSummaries,
804 instrumentName="",
805 refEpoch=None,
806 refObjectLoader=None,
807 inputCameraModel=None,
808 colorCatalog=None,
809 inputCamera=None,
810 nCores=1,
811 ):
812 """Run the WCS fit for a given set of visits
814 Parameters
815 ----------
816 inputCatalogRefs : `list` [`DeferredDatasetHandle`]
817 List of handles pointing to visit-level source
818 tables.
819 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
820 List of catalogs with per-detector summary information.
821 instrumentName : `str`, optional
822 Name of the instrument used. This is only used for labelling.
823 refEpoch : `float`
824 Epoch of the reference objects in MJD.
825 refObjectLoader : instance of
826 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`
827 Reference object loader instance.
828 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
829 Parameters to use for the device part of the model.
830 colorCatalog : `lsst.afw.table.SimpleCatalog`
831 Catalog containing object coordinates and magnitudes.
832 inputCamera : `lsst.afw.cameraGeom.Camera`, optional
833 Camera to be used as template when constructing new camera.
834 nCores : `int`, optional
835 Number of cores to use during WCS fit.
837 Returns
838 -------
839 result : `lsst.pipe.base.Struct`
840 ``outputWcss`` : `list` [`lsst.afw.table.ExposureCatalog`]
841 List of exposure catalogs (one per visit) with the WCS for each
842 detector set by the new fitted WCS.
843 ``fitModel`` : `wcsfit.WCSFit`
844 Model-fitting object with final model parameters.
845 ``outputCatalog`` : `pyarrow.Table`
846 Catalog with fit residuals of all sources used.
847 ``starCatalog`` : `pyarrow.Table`
848 Catalog with best-fit positions of the objects fit.
849 ``modelParams`` : `dict`
850 Parameters and covariance of the best-fit WCS model.
851 ``cameraModelParams`` : `dict` [`str`, `np.ndarray`]
852 Parameters of the device part of the model, in the format
853 needed as input for future runs.
854 ``colorParams`` : `dict` [`int`, `np.ndarray`]
855 DCR parameters fit in RA and Dec directions for each visit.
856 ``camera`` : `lsst.afw.cameraGeom.Camera`
857 Camera object constructed from the per-detector model.
858 """
859 self.log.info("Gather instrument, exposure, and field info")
861 # Get RA, Dec, MJD, etc., for the input visits
862 exposureInfo, exposuresHelper, extensionInfo, instruments = self._get_exposure_info(
863 inputVisitSummaries
864 )
866 # Get information about the extent of the input visits
867 fields, fieldCenter, fieldRadius = self._prep_sky(inputVisitSummaries, exposureInfo.medianEpoch)
868 self.log.info("Field center set at %s with radius %s degrees", fieldCenter, fieldRadius.asDegrees())
870 self.log.info("Load catalogs and associate sources")
871 # Set up class to associate sources into matches using a
872 # friends-of-friends algorithm
873 associations = wcsfit.FoFClass(
874 fields,
875 instruments,
876 exposuresHelper,
877 [fieldRadius.asDegrees()],
878 (self.config.matchRadius * u.arcsec).to(u.degree).value,
879 )
881 # Add the reference catalog to the associator
882 medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="jyear").mjd
883 refObjects, refCovariance = self._load_refcat(
884 refObjectLoader,
885 extensionInfo,
886 epoch=medianEpoch,
887 center=fieldCenter,
888 radius=fieldRadius,
889 associations=associations,
890 )
892 # Add the science catalogs and associate new sources as they are added
893 sourceIndices, usedColumns = self._load_catalogs_and_associate(
894 associations, inputCatalogRefs, extensionInfo
895 )
896 self._check_degeneracies(associations, extensionInfo)
898 self.log.info("Fit the WCSs")
899 # Set up a YAML-type string using the config variables and a sample
900 # visit
901 inputYaml, mapTemplate = self.make_yaml(
902 inputVisitSummaries[0],
903 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
904 )
906 # Set the verbosity level for WCSFit from the task log level.
907 # TODO: DM-36850, Add lsst.log to gbdes so that log messages are
908 # properly propagated.
909 loglevel = self.log.getEffectiveLevel()
910 if loglevel >= self.log.WARNING:
911 verbose = 0
912 elif loglevel == self.log.INFO:
913 verbose = 1
914 else:
915 verbose = 2
917 # Set up the WCS-fitting class using the results of the FOF associator
918 if self.config.useInputCameraModel:
919 fixMaps = ",".join(inputCameraModel.keys())
920 else:
921 fixMaps = ""
922 wcsf = wcsfit.WCSFit(
923 fields,
924 instruments,
925 exposuresHelper,
926 extensionInfo.visitIndex,
927 extensionInfo.detectorIndex,
928 inputYaml,
929 extensionInfo.wcs,
930 associations.sequence,
931 associations.extn,
932 associations.obj,
933 sysErr=self.config.systematicError,
934 refSysErr=self.config.referenceSystematicError,
935 usePM=self.config.fitProperMotion,
936 verbose=verbose,
937 fixMaps=fixMaps,
938 num_threads=nCores,
939 )
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").jyear
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 neededColumns = ["coord_ra", "coord_dec", "coord_raErr", "coord_decErr"]
1329 if self.config.applyRefCatProperMotion:
1330 neededColumns += [
1331 "pm_ra",
1332 "pm_dec",
1333 "parallax",
1334 "pm_raErr",
1335 "pm_decErr",
1336 "parallaxErr",
1337 ]
1338 # Get refcat version from refcat metadata
1339 refCatMetadata = refObjectLoader.refCats[0].get().getMetadata()
1340 # DM-47181: Added this to work for LSSTComCam with
1341 # the_monster_20240904 catalog that does not have this key.
1342 if "REFCAT_FORMAT_VERSION" not in refCatMetadata:
1343 refCatVersion = 2
1344 else:
1345 refCatVersion = refCatMetadata["REFCAT_FORMAT_VERSION"]
1346 if refCatVersion == 2:
1347 neededColumns += [
1348 "coord_ra_coord_dec_Cov",
1349 "coord_ra_pm_ra_Cov",
1350 "coord_ra_pm_dec_Cov",
1351 "coord_ra_parallax_Cov",
1352 "coord_dec_pm_ra_Cov",
1353 "coord_dec_pm_dec_Cov",
1354 "coord_dec_parallax_Cov",
1355 "pm_ra_pm_dec_Cov",
1356 "pm_ra_parallax_Cov",
1357 "pm_dec_parallax_Cov",
1358 ]
1360 # Load each shard of the reference catalog separately to avoid a spike
1361 # in the memory load.
1362 refCatShards = []
1363 for dataId, cat in zip(refObjectLoader.dataIds, refObjectLoader.refCats):
1364 miniRefObjectLoader = ReferenceObjectLoader(
1365 dataIds=[dataId],
1366 refCats=[cat],
1367 config=refObjectLoader.config,
1368 log=self.log,
1369 )
1370 try:
1371 if region is not None:
1372 skyRegion = miniRefObjectLoader.loadRegion(
1373 region, self.config.referenceFilter, epoch=formattedEpoch
1374 )
1375 elif (center is not None) and (radius is not None):
1376 skyRegion = miniRefObjectLoader.loadSkyCircle(
1377 center, radius, self.config.referenceFilter, epoch=formattedEpoch
1378 )
1379 else:
1380 raise RuntimeError("Either `region` or `center` and `radius` must be set.")
1381 except RuntimeError:
1382 self.log.debug("Reference catalog shard has no objects inside the region.")
1383 continue
1384 selected = self.referenceSelector.run(skyRegion.refCat)
1385 # Need memory contiguity to get reference filters as a vector.
1386 if not selected.sourceCat.isContiguous():
1387 refCatShard = selected.sourceCat.copy(deep=True)
1388 else:
1389 refCatShard = selected.sourceCat
1390 refCatShard = refCatShard.asAstropy()[neededColumns]
1392 # In Gaia DR3, missing values are denoted by NaNs.
1393 finiteInd = np.isfinite(refCatShard["coord_ra"]) & np.isfinite(refCatShard["coord_dec"])
1394 refCatShard = refCatShard[finiteInd]
1395 refCatShards.append(refCatShard)
1396 refCat = vstack(refCatShards)
1398 if self.config.excludeNonPMObjects and self.config.applyRefCatProperMotion:
1399 # Gaia DR2 has zeros for missing data, while Gaia DR3 has NaNs:
1400 hasPM = (
1401 (refCat["pm_raErr"] != 0) & np.isfinite(refCat["pm_raErr"]) & np.isfinite(refCat["pm_decErr"])
1402 )
1403 refCat = refCat[hasPM]
1405 ra = (refCat["coord_ra"]).to(u.degree).to_value().tolist()
1406 dec = (refCat["coord_dec"]).to(u.degree).to_value().tolist()
1407 raCov = ((refCat["coord_raErr"]).to(u.degree).to_value() ** 2).tolist()
1408 decCov = ((refCat["coord_decErr"]).to(u.degree).to_value() ** 2).tolist()
1410 if refCatVersion == 2:
1411 raDecCov = (refCat["coord_ra_coord_dec_Cov"]).to(u.degree**2).to_value().tolist()
1412 else:
1413 raDecCov = np.zeros(len(ra))
1415 refObjects = {"ra": ra, "dec": dec, "raCov": raCov, "decCov": decCov, "raDecCov": raDecCov}
1416 refCovariance = []
1418 if self.config.fitProperMotion:
1419 raPM = (refCat["pm_ra"]).to(u.marcsec).to_value().tolist()
1420 decPM = (refCat["pm_dec"]).to(u.marcsec).to_value().tolist()
1421 parallax = (refCat["parallax"]).to(u.marcsec).to_value().tolist()
1422 cov = _make_ref_covariance_matrix(refCat, version=refCatVersion)
1423 pmDict = {"raPM": raPM, "decPM": decPM, "parallax": parallax}
1424 refObjects.update(pmDict)
1425 refCovariance = cov
1427 if associations is not None:
1428 extensionIndex = np.flatnonzero(extensionInfo.extensionType == "REFERENCE")[0]
1429 visitIndex = extensionInfo.visitIndex[extensionIndex]
1430 detectorIndex = extensionInfo.detectorIndex[extensionIndex]
1431 instrumentIndex = -1 # -1 indicates the reference catalog
1432 refWcs = extensionInfo.wcs[extensionIndex]
1434 associations.addCatalog(
1435 refWcs,
1436 "STELLAR",
1437 visitIndex,
1438 fieldIndex,
1439 instrumentIndex,
1440 detectorIndex,
1441 extensionIndex,
1442 np.ones(len(refCat), dtype=bool),
1443 ra,
1444 dec,
1445 np.arange(len(ra)),
1446 )
1448 return refObjects, refCovariance
1450 @staticmethod
1451 def _find_extension_index(extensionInfo, visit, detector):
1452 """Find the index for a given extension from its visit and detector
1453 number.
1455 If no match is found, None is returned.
1457 Parameters
1458 ----------
1459 extensionInfo : `lsst.pipe.base.Struct`
1460 Struct containing properties for each extension.
1461 visit : `int`
1462 Visit number
1463 detector : `int`
1464 Detector number
1466 Returns
1467 -------
1468 extensionIndex : `int` or None
1469 Index of this extension
1470 """
1471 findExtension = np.flatnonzero((extensionInfo.visit == visit) & (extensionInfo.detector == detector))
1472 if len(findExtension) == 0:
1473 extensionIndex = None
1474 else:
1475 extensionIndex = findExtension[0]
1476 return extensionIndex
1478 def _load_catalogs_and_associate(
1479 self, associations, inputCatalogRefs, extensionInfo, fieldIndex=0, instrumentIndex=0
1480 ):
1481 """Load the science catalogs and add the sources to the associator
1482 class `wcsfit.FoFClass`, associating them into matches as you go.
1484 Parameters
1485 ----------
1486 associations : `wcsfit.FoFClass`
1487 Object to which to add the catalog of source and which performs
1488 the source association.
1489 inputCatalogRefs : `list`
1490 List of DeferredDatasetHandles pointing to visit-level source
1491 tables.
1492 extensionInfo : `lsst.pipe.base.Struct`
1493 Struct containing properties for each extension (visit/detector).
1494 ``visit`` : `np.ndarray`
1495 Name of the visit for this extension.
1496 ``detector`` : `np.ndarray`
1497 Name of the detector for this extension.
1498 ``visitIndex` : `np.ndarray` [`int`]
1499 Index of visit for this extension.
1500 ``detectorIndex`` : `np.ndarray` [`int`]
1501 Index of the detector for this extension.
1502 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1503 Initial WCS for this extension.
1504 ``extensionType`` : `np.ndarray` [`str`]
1505 "SCIENCE" or "REFERENCE".
1506 fieldIndex : `int`
1507 Index of the field for these catalogs. Should be zero assuming all
1508 data is being fit together.
1509 instrumentIndex : `int`
1510 Index of the instrument for these catalogs. Should be zero
1511 assuming all data comes from the same instrument.
1513 Returns
1514 -------
1515 sourceIndices : `list`
1516 List of boolean arrays used to select sources.
1517 columns : `list` [`str`]
1518 List of columns needed from source tables.
1519 """
1520 columns = [
1521 "detector",
1522 "sourceId",
1523 "x",
1524 "xErr",
1525 "y",
1526 "yErr",
1527 "ixx",
1528 "iyy",
1529 "ixy",
1530 f"{self.config.sourceFluxType}_instFlux",
1531 f"{self.config.sourceFluxType}_instFluxErr",
1532 ]
1533 if self.sourceSelector.config.doFlags:
1534 columns.extend(self.sourceSelector.config.flags.bad)
1535 if self.sourceSelector.config.doUnresolved:
1536 columns.append(self.sourceSelector.config.unresolved.name)
1537 if self.sourceSelector.config.doIsolated:
1538 columns.append(self.sourceSelector.config.isolated.parentName)
1539 columns.append(self.sourceSelector.config.isolated.nChildName)
1540 if self.sourceSelector.config.doRequirePrimary:
1541 columns.append(self.sourceSelector.config.requirePrimary.primaryColName)
1543 sourceIndices = [None] * len(extensionInfo.visit)
1544 for inputCatalogRef in inputCatalogRefs:
1545 visit = inputCatalogRef.dataId["visit"]
1546 inputCatalog = inputCatalogRef.get(parameters={"columns": columns})
1547 # Get a sorted array of detector names
1548 detectors = np.unique(inputCatalog["detector"])
1550 for detector in detectors:
1551 detectorSources = inputCatalog[inputCatalog["detector"] == detector]
1552 xCov = detectorSources["xErr"] ** 2
1553 yCov = detectorSources["yErr"] ** 2
1554 xyCov = (
1555 detectorSources["ixy"] * (xCov + yCov) / (detectorSources["ixx"] + detectorSources["iyy"])
1556 )
1557 # Remove sources with bad shape measurements
1558 goodShapes = xyCov**2 <= (xCov * yCov)
1559 selected = self.sourceSelector.run(detectorSources)
1560 goodInds = selected.selected & goodShapes
1561 isStar = np.ones(goodInds.sum())
1562 extensionIndex = self._find_extension_index(extensionInfo, visit, detector)
1563 if extensionIndex is None:
1564 # This extension does not have information necessary for
1565 # fit. Skip the detections from this detector for this
1566 # visit.
1567 continue
1568 detectorIndex = extensionInfo.detectorIndex[extensionIndex]
1569 visitIndex = extensionInfo.visitIndex[extensionIndex]
1571 sourceIndices[extensionIndex] = goodInds
1573 wcs = extensionInfo.wcs[extensionIndex]
1574 associations.reprojectWCS(wcs, fieldIndex)
1576 associations.addCatalog(
1577 wcs,
1578 "STELLAR",
1579 visitIndex,
1580 fieldIndex,
1581 instrumentIndex,
1582 detectorIndex,
1583 extensionIndex,
1584 isStar,
1585 detectorSources[goodInds]["x"].to_list(),
1586 detectorSources[goodInds]["y"].to_list(),
1587 np.arange(goodInds.sum()),
1588 )
1590 associations.sortMatches(
1591 fieldIndex, minMatches=self.config.minMatches, allowSelfMatches=self.config.allowSelfMatches
1592 )
1594 return sourceIndices, columns
1596 def _check_degeneracies(self, associations, extensionInfo):
1597 """Check that the minimum number of visits and sources needed to
1598 constrain the model are present.
1600 This does not guarantee that the Hessian matrix of the chi-square,
1601 which is used to fit the model, will be positive-definite, but if the
1602 checks here do not pass, the matrix is certain to not be
1603 positive-definite and the model cannot be fit.
1605 Parameters
1606 ----------
1607 associations : `wcsfit.FoFClass`
1608 Object holding the source association information.
1609 extensionInfo : `lsst.pipe.base.Struct`
1610 Struct containing properties for each extension (visit/detector):
1611 ``visit`` : `np.ndarray`
1612 Name of the visit for this extension.
1613 ``detector`` : `np.ndarray`
1614 Name of the detector for this extension.
1615 ``visitIndex` : `np.ndarray` [`int`]
1616 Index of visit for this extension.
1617 ``detectorIndex`` : `np.ndarray` [`int`]
1618 Index of the detector for this extension.
1619 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1620 Initial WCS for this extension.
1621 ``extensionType`` : `np.ndarray` [`str`]
1622 "SCIENCE" or "REFERENCE".
1623 """
1624 # As a baseline, need to have more stars per detector than per-detector
1625 # parameters, and more stars per visit than per-visit parameters.
1626 whichExtension = np.array(associations.extn)
1627 whichDetector = np.zeros(len(whichExtension))
1628 whichVisit = np.zeros(len(whichExtension))
1630 for extension, (detector, visit) in enumerate(zip(extensionInfo.detector, extensionInfo.visit)):
1631 ex_ind = whichExtension == extension
1632 whichDetector[ex_ind] = detector
1633 whichVisit[ex_ind] = visit
1635 if (not self.config.useInputCameraModel) and ("BAND/DEVICE/poly" in self.config.deviceModel):
1636 nCoeffDetectorModel = _nCoeffsFromDegree(self.config.devicePolyOrder)
1637 unconstrainedDetectors = []
1638 for detector in np.unique(extensionInfo.detector):
1639 numSources = (whichDetector == detector).sum()
1640 if numSources < nCoeffDetectorModel:
1641 unconstrainedDetectors.append(str(detector))
1643 if unconstrainedDetectors:
1644 raise pipeBase.NoWorkFound(
1645 "The model is not constrained. The following detectors do not have enough "
1646 f"sources ({nCoeffDetectorModel} required): ",
1647 ", ".join(unconstrainedDetectors),
1648 )
1650 def make_yaml(self, inputVisitSummary, inputFile=None, inputCameraModel=None):
1651 """Make a YAML-type object that describes the parameters of the fit
1652 model.
1654 Parameters
1655 ----------
1656 inputVisitSummary : `lsst.afw.table.ExposureCatalog`
1657 Catalog with per-detector summary information.
1658 inputFile : `str`
1659 Path to a file that contains a basic model.
1660 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
1661 Parameters to use for the device part of the model.
1663 Returns
1664 -------
1665 inputYaml : `wcsfit.YAMLCollector`
1666 YAML object containing the model description.
1667 inputDict : `dict` [`str`, `str`]
1668 Dictionary containing the model description.
1669 """
1670 if inputFile is not None:
1671 inputYaml = wcsfit.YAMLCollector(inputFile, "PixelMapCollection")
1672 else:
1673 inputYaml = wcsfit.YAMLCollector("", "PixelMapCollection")
1674 inputDict = {}
1675 modelComponents = ["BAND/DEVICE", "EXPOSURE"]
1676 baseMap = {"Type": "Composite", "Elements": modelComponents}
1677 inputDict["EXPOSURE/DEVICE/base"] = baseMap
1679 xMin = str(inputVisitSummary["bbox_min_x"].min())
1680 xMax = str(inputVisitSummary["bbox_max_x"].max())
1681 yMin = str(inputVisitSummary["bbox_min_y"].min())
1682 yMax = str(inputVisitSummary["bbox_max_y"].max())
1684 deviceModel = {"Type": "Composite", "Elements": self.config.deviceModel.list()}
1685 inputDict["BAND/DEVICE"] = deviceModel
1686 for component in self.config.deviceModel:
1687 mapping = component.split("/")[-1]
1688 if mapping == "poly":
1689 componentDict = {
1690 "Type": "Poly",
1691 "XPoly": {"OrderX": self.config.devicePolyOrder, "SumOrder": True},
1692 "YPoly": {"OrderX": self.config.devicePolyOrder, "SumOrder": True},
1693 "XMin": xMin,
1694 "XMax": xMax,
1695 "YMin": yMin,
1696 "YMax": yMax,
1697 }
1698 elif mapping == "identity":
1699 componentDict = {"Type": "Identity"}
1701 inputDict[component] = componentDict
1703 if (inputCameraModel is not None) and self.config.useInputCameraModel:
1704 # This assumes that the input camera model is a 'poly' model
1705 nCoeffs = _nCoeffsFromDegree(self.config.devicePolyOrder)
1706 for key, coeffs in inputCameraModel.items():
1707 if len(coeffs) != nCoeffs * 2:
1708 raise RuntimeError(
1709 "Input camera model polynomial order does not match the devicePolyOrder"
1710 )
1711 mapDict = {
1712 "Type": "Poly",
1713 "XPoly": {
1714 "OrderX": self.config.devicePolyOrder,
1715 "SumOrder": True,
1716 "Coefficients": coeffs[:nCoeffs].tolist(),
1717 },
1718 "YPoly": {
1719 "OrderX": self.config.devicePolyOrder,
1720 "SumOrder": True,
1721 "Coefficients": coeffs[nCoeffs:].tolist(),
1722 },
1723 "XMin": xMin,
1724 "XMax": xMax,
1725 "YMin": yMin,
1726 "YMax": yMax,
1727 }
1728 inputDict[key] = mapDict
1730 exposureModelComponents = self.config.exposureModel.list()
1731 if self.config.useColor:
1732 exposureModelComponents.append("EXPOSURE/dcr")
1733 exposureModel = {"Type": "Composite", "Elements": exposureModelComponents}
1734 inputDict["EXPOSURE"] = exposureModel
1735 for component in exposureModelComponents:
1736 mapping = component.split("/")[-1]
1737 if mapping == "poly":
1738 componentDict = {
1739 "Type": "Poly",
1740 "XPoly": {"OrderX": self.config.exposurePolyOrder, "SumOrder": "true"},
1741 "YPoly": {"OrderX": self.config.exposurePolyOrder, "SumOrder": "true"},
1742 }
1743 elif mapping == "identity":
1744 componentDict = {"Type": "Identity"}
1745 elif mapping == "dcr":
1746 componentDict = {
1747 "Type": "Color",
1748 "Reference": self.config.referenceColor,
1749 "Function": {"Type": "Constant"},
1750 }
1752 inputDict[component] = componentDict
1754 inputYaml.addInput(yaml.dump(inputDict))
1755 inputYaml.addInput("Identity:\n Type: Identity\n")
1757 return inputYaml, inputDict
1759 def _add_objects(self, wcsf, inputCatalogRefs, sourceIndices, extensionInfo, columns):
1760 """Add science sources to the wcsfit.WCSFit object.
1762 Parameters
1763 ----------
1764 wcsf : `wcsfit.WCSFit`
1765 WCS-fitting object.
1766 inputCatalogRefs : `list`
1767 List of DeferredDatasetHandles pointing to visit-level source
1768 tables.
1769 sourceIndices : `list`
1770 List of boolean arrays used to select sources.
1771 extensionInfo : `lsst.pipe.base.Struct`
1772 Struct containing properties for each extension (visit/detector):
1773 ``visit`` : `np.ndarray`
1774 Name of the visit for this extension.
1775 ``detector`` : `np.ndarray`
1776 Name of the detector for this extension.
1777 ``visitIndex` : `np.ndarray` [`int`]
1778 Index of visit for this extension.
1779 ``detectorIndex`` : `np.ndarray` [`int`]
1780 Index of the detector for this extension.
1781 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1782 Initial WCS for this extension.
1783 ``extensionType`` : `np.ndarray` [`str`]
1784 "SCIENCE" or "REFERENCE".
1785 columns : `list` [`str`]
1786 List of columns needed from source tables.
1787 """
1788 for inputCatalogRef in inputCatalogRefs:
1789 visit = inputCatalogRef.dataId["visit"]
1790 inputCatalog = inputCatalogRef.get(parameters={"columns": columns})
1791 detectors = np.unique(inputCatalog["detector"])
1793 for detector in detectors:
1794 detectorSources = inputCatalog[inputCatalog["detector"] == detector]
1796 extensionIndex = self._find_extension_index(extensionInfo, visit, detector)
1797 if extensionIndex is None:
1798 # This extension does not have information necessary for
1799 # fit. Skip the detections from this detector for this
1800 # visit.
1801 continue
1803 sourceCat = detectorSources[sourceIndices[extensionIndex]]
1805 xCov = sourceCat["xErr"] ** 2
1806 yCov = sourceCat["yErr"] ** 2
1807 xyCov = sourceCat["ixy"] * (xCov + yCov) / (sourceCat["ixx"] + sourceCat["iyy"])
1808 # TODO: add correct xyErr if DM-7101 is ever done.
1810 d = {
1811 "x": sourceCat["x"].to_numpy(),
1812 "y": sourceCat["y"].to_numpy(),
1813 "xCov": xCov.to_numpy(),
1814 "yCov": yCov.to_numpy(),
1815 "xyCov": xyCov.to_numpy(),
1816 }
1818 wcsf.setObjects(
1819 extensionIndex,
1820 d,
1821 "x",
1822 "y",
1823 ["xCov", "yCov", "xyCov"],
1824 defaultColor=self.config.referenceColor,
1825 )
1827 def _add_ref_objects(self, wcsf, refObjects, refCovariance, extensionInfo, fieldIndex=0):
1828 """Add reference sources to the wcsfit.WCSFit object.
1830 Parameters
1831 ----------
1832 wcsf : `wcsfit.WCSFit`
1833 WCS-fitting object.
1834 refObjects : `dict`
1835 Position and error information of reference objects.
1836 refCovariance : `list` [`float`]
1837 Flattened output covariance matrix.
1838 extensionInfo : `lsst.pipe.base.Struct`
1839 Struct containing properties for each extension (visit/detector):
1840 ``visit`` : `np.ndarray`
1841 Name of the visit for this extension.
1842 ``detector`` : `np.ndarray`
1843 Name of the detector for this extension.
1844 ``visitIndex` : `np.ndarray` [`int`]
1845 Index of visit for this extension.
1846 ``detectorIndex`` : `np.ndarray` [`int`]
1847 Index of the detector for this extension.
1848 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
1849 Initial WCS for this extension.
1850 ``extensionType`` : `np.ndarray` [`str`]
1851 "SCIENCE" or "REFERENCE".
1852 fieldIndex : `int`, optional
1853 Index of the field to which these sources correspond.
1854 """
1855 extensionIndex = np.flatnonzero(
1856 (extensionInfo.extensionType == "REFERENCE") & (extensionInfo.visit == fieldIndex)
1857 )[0]
1858 if self.config.fitProperMotion:
1859 wcsf.setObjects(
1860 extensionIndex,
1861 refObjects,
1862 "ra",
1863 "dec",
1864 ["raCov", "decCov", "raDecCov"],
1865 pmDecKey="decPM",
1866 pmRaKey="raPM",
1867 parallaxKey="parallax",
1868 pmCovKey="fullCov",
1869 pmCov=refCovariance,
1870 )
1871 else:
1872 wcsf.setObjects(extensionIndex, refObjects, "ra", "dec", ["raCov", "decCov", "raDecCov"])
1874 def _add_color_objects(self, wcsf, colorCatalog):
1875 """Associate input matches with objects in color catalog and set their
1876 color value.
1878 Parameters
1879 ----------
1880 wcsf : `wcsfit.WCSFit`
1881 WCSFit object, assumed to have fit model.
1882 colorCatalog : `lsst.afw.table.SimpleCatalog`
1883 Catalog containing object coordinates and magnitudes.
1884 """
1886 # Get current best position for matches
1887 starCat = wcsf.getStarCatalog()
1889 # TODO: DM-45650, update how the colors are read in here.
1890 catalogBands = colorCatalog.metadata.getArray("BANDS")
1891 colorInd1 = catalogBands.index(self.config.color[0])
1892 colorInd2 = catalogBands.index(self.config.color[1])
1893 colors = colorCatalog["mag_std_noabs"][:, colorInd1] - colorCatalog["mag_std_noabs"][:, colorInd2]
1894 goodInd = (colorCatalog["mag_std_noabs"][:, colorInd1] != 99.0) & (
1895 colorCatalog["mag_std_noabs"][:, colorInd2] != 99.0
1896 )
1898 with Matcher(np.array(starCat["starX"]), np.array(starCat["starY"])) as matcher:
1899 idx, idx_starCat, idx_colorCat, d = matcher.query_radius(
1900 (colorCatalog[goodInd]["coord_ra"] * u.radian).to(u.degree).value,
1901 (colorCatalog[goodInd]["coord_dec"] * u.radian).to(u.degree).value,
1902 self.config.matchRadius / 3600.0,
1903 return_indices=True,
1904 )
1906 matchesWithColor = starCat["starMatchID"][idx_starCat]
1907 matchColors = np.ones(len(matchesWithColor)) * self.config.referenceColor
1908 matchColors = colors[goodInd][idx_colorCat]
1909 wcsf.setColors(matchesWithColor, matchColors)
1911 def _make_afw_wcs(self, mapDict, centerRA, centerDec, doNormalizePixels=False, xScale=1, yScale=1):
1912 """Make an `lsst.afw.geom.SkyWcs` from a dictionary of mappings.
1914 Parameters
1915 ----------
1916 mapDict : `dict`
1917 Dictionary of mapping parameters.
1918 centerRA : `lsst.geom.Angle`
1919 RA of the tangent point.
1920 centerDec : `lsst.geom.Angle`
1921 Declination of the tangent point.
1922 doNormalizePixels : `bool`
1923 Whether to normalize pixels so that range is [-1,1].
1924 xScale : `float`
1925 Factor by which to normalize x-dimension. Corresponds to width of
1926 detector.
1927 yScale : `float`
1928 Factor by which to normalize y-dimension. Corresponds to height of
1929 detector.
1931 Returns
1932 -------
1933 outWCS : `lsst.afw.geom.SkyWcs`
1934 WCS constructed from the input mappings
1935 """
1936 # Set up pixel frames
1937 pixelFrame = astshim.Frame(2, "Domain=PIXELS")
1938 normedPixelFrame = astshim.Frame(2, "Domain=NORMEDPIXELS")
1940 if doNormalizePixels:
1941 # Pixels will need to be rescaled before going into the mappings
1942 normCoefficients = [-1.0, 2.0 / xScale, 0, -1.0, 0, 2.0 / yScale]
1943 normMap = _convert_to_ast_polymap_coefficients(normCoefficients)
1944 else:
1945 normMap = astshim.UnitMap(2)
1947 # All of the detectors for one visit map to the same tangent plane
1948 tangentPoint = lsst.geom.SpherePoint(centerRA, centerDec)
1949 cdMatrix = afwgeom.makeCdMatrix(1.0 * lsst.geom.degrees, 0 * lsst.geom.degrees, True)
1950 iwcToSkyWcs = afwgeom.makeSkyWcs(lsst.geom.Point2D(0, 0), tangentPoint, cdMatrix)
1951 iwcToSkyMap = iwcToSkyWcs.getFrameDict().getMapping("PIXELS", "SKY")
1952 skyFrame = iwcToSkyWcs.getFrameDict().getFrame("SKY")
1954 frameDict = astshim.FrameDict(pixelFrame)
1955 frameDict.addFrame("PIXELS", normMap, normedPixelFrame)
1957 currentFrameName = "NORMEDPIXELS"
1959 # Dictionary values are ordered according to the maps' application.
1960 for m, mapElement in enumerate(mapDict.values()):
1961 mapType = mapElement["Type"]
1963 if mapType == "Poly":
1964 mapCoefficients = mapElement["Coefficients"]
1965 astMap = _convert_to_ast_polymap_coefficients(mapCoefficients)
1966 elif mapType == "Identity":
1967 astMap = astshim.UnitMap(2)
1968 else:
1969 raise ValueError(f"Converting map type {mapType} to WCS is not supported")
1971 if m == len(mapDict) - 1:
1972 newFrameName = "IWC"
1973 else:
1974 newFrameName = "INTERMEDIATE" + str(m)
1975 newFrame = astshim.Frame(2, f"Domain={newFrameName}")
1976 frameDict.addFrame(currentFrameName, astMap, newFrame)
1977 currentFrameName = newFrameName
1978 frameDict.addFrame("IWC", iwcToSkyMap, skyFrame)
1980 outWCS = afwgeom.SkyWcs(frameDict)
1981 return outWCS
1983 def _make_outputs(
1984 self,
1985 wcsf,
1986 visitSummaryTables,
1987 exposureInfo,
1988 mapTemplate,
1989 extensionInfo,
1990 inputCameraModel=None,
1991 inputCamera=None,
1992 ):
1993 """Make a WCS object out of the WCS models.
1995 Parameters
1996 ----------
1997 wcsf : `wcsfit.WCSFit`
1998 WCSFit object, assumed to have fit model.
1999 visitSummaryTables : `list` [`lsst.afw.table.ExposureCatalog`]
2000 Catalogs with per-detector summary information from which to grab
2001 detector information.
2002 extensionInfo : `lsst.pipe.base.Struct`
2003 Struct containing properties for each extension (visit/detector):
2004 ``visit`` : `np.ndarray`
2005 Name of the visit for this extension.
2006 ``detector`` : `np.ndarray`
2007 Name of the detector for this extension.
2008 ``visitIndex` : `np.ndarray` [`int`]
2009 Index of visit for this extension.
2010 ``detectorIndex`` : `np.ndarray` [`int`]
2011 Index of the detector for this extension.
2012 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2013 Initial WCS for this extension.
2014 ``extensionType`` : `np.ndarray` [`str`]
2015 "SCIENCE" or "REFERENCE".
2016 mapTemplate : `dict` [`str`, `str`]
2017 Dictionary containing the model description.
2018 extensionInfo : `lsst.pipe.base.Struct`
2019 Struct containing properties for each extension (visit/detector).
2020 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
2021 Parameters to use for the device part of the model. This must be
2022 provided if an input camera model was used.
2024 Returns
2025 -------
2026 catalogs : `dict` [`str`, `lsst.afw.table.ExposureCatalog`]
2027 Dictionary of `lsst.afw.table.ExposureCatalog` objects with the WCS
2028 set to the WCS fit in wcsf, keyed by the visit name.
2029 cameraParams : `dict` [`str`, `np.ndarray`], optional
2030 Parameters for the device part of the model in the format needed
2031 when used as input for future runs.
2032 colorFits : `dict` [`int`, `np.ndarray`], optional
2033 DCR parameters fit in RA and Dec directions for each visit.
2034 camera : `lsst.afw.cameraGeom.Camera`, optional
2035 Camera object constructed from the per-detector model.
2036 """
2037 # Get the parameters of the fit models
2038 mapParams = wcsf.mapCollection.getParamDict()
2039 cameraParams = {}
2040 if self.config.saveCameraModel:
2041 for element in mapTemplate["BAND/DEVICE"]["Elements"]:
2042 for detector in exposureInfo.detectors:
2043 detectorTemplate = element.replace("DEVICE", str(detector))
2044 detectorTemplate = detectorTemplate.replace("BAND", ".+")
2045 for k, params in mapParams.items():
2046 if re.fullmatch(detectorTemplate, k):
2047 cameraParams[k] = params
2048 if self.config.saveCameraObject:
2049 # Get the average rotation angle of the input visits.
2050 rotations = [
2051 visTable[0].visitInfo.boresightRotAngle.asRadians() for visTable in visitSummaryTables
2052 ]
2053 rotationAngle = np.mean(rotations)
2054 if inputCamera is None:
2055 raise RuntimeError(
2056 "inputCamera must be provided to _make_outputs in order to build output camera."
2057 )
2058 camera = self.buildCamera.run(
2059 mapParams,
2060 mapTemplate,
2061 exposureInfo.detectors,
2062 exposureInfo.visits,
2063 inputCamera,
2064 rotationAngle,
2065 )
2066 else:
2067 camera = None
2068 if self.config.useInputCameraModel:
2069 if inputCameraModel is None:
2070 raise RuntimeError(
2071 "inputCameraModel must be provided to _make_outputs in order to build output WCS."
2072 )
2073 mapParams.update(inputCameraModel)
2075 # Set up the schema for the output catalogs
2076 schema = lsst.afw.table.ExposureTable.makeMinimalSchema()
2077 schema.addField("visit", type="L", doc="Visit number")
2078 schema.addField(
2079 "recoveredWcs",
2080 type="Flag",
2081 doc="Input WCS missing, output recovered from other input visit/detectors.",
2082 )
2084 # Pixels will need to be rescaled before going into the mappings
2085 xscale = int(mapTemplate["BAND/DEVICE/poly"]["XMax"]) - int(mapTemplate["BAND/DEVICE/poly"]["XMin"])
2086 yscale = int(mapTemplate["BAND/DEVICE/poly"]["YMax"]) - int(mapTemplate["BAND/DEVICE/poly"]["YMin"])
2088 catalogs = {}
2089 colorFits = {}
2090 partialOutputs = False
2091 for v, visitSummary in enumerate(visitSummaryTables):
2092 visit = visitSummary[0]["visit"]
2093 if visit not in extensionInfo.visit:
2094 self.log.warning("Visit %d was dropped because no detectors had valid WCSs.", visit)
2095 partialOutputs = True
2096 continue
2098 visitMaps = wcsf.mapCollection.orderAtoms(f"{visit}")
2099 if self.config.useColor:
2100 colorMap = visitMaps.pop(visitMaps.index(f"{visit}/dcr"))
2101 colorFits[visit] = mapParams[colorMap]
2102 visitMap = visitMaps[0]
2103 visitMapType = wcsf.mapCollection.getMapType(visitMap)
2104 if (visitMap not in mapParams) and (visitMapType != "Identity"):
2105 self.log.warning("Visit %d was dropped because of an insufficient amount of data.", visit)
2106 partialOutputs = True
2107 continue
2109 catalog = lsst.afw.table.ExposureCatalog(schema)
2110 catalog.resize(len(exposureInfo.detectors))
2111 catalog["visit"] = visit
2113 for d, detector in enumerate(exposureInfo.detectors):
2114 mapName = f"{visit}/{detector}"
2115 if mapName in wcsf.mapCollection.allMapNames():
2116 mapElements = wcsf.mapCollection.orderAtoms(f"{mapName}/base")
2117 catalog[d]["recoveredWcs"] = False
2118 else:
2119 # This extension was not fit, but the WCS can be recovered
2120 # using the maps fit from sources on other visits but the
2121 # same detector and from sources on other detectors from
2122 # this visit.
2123 genericElements = mapTemplate["EXPOSURE/DEVICE/base"]["Elements"]
2124 mapElements = []
2125 band = visitSummary[0]["physical_filter"]
2127 # Go through the generic map components to build the names
2128 # of the specific maps for this extension.
2129 for component in genericElements:
2130 elements = mapTemplate[component]["Elements"]
2131 for element in elements:
2132 # TODO: DM-42519, gbdes sets the "BAND" to the
2133 # instrument name currently. This will need to be
2134 # disambiguated if we run on multiple bands at
2135 # once.
2136 element = element.replace("BAND", str(band))
2137 element = element.replace("EXPOSURE", str(visit))
2138 element = element.replace("DEVICE", str(detector))
2139 mapElements.append(element)
2140 catalog[d]["recoveredWcs"] = True
2141 mapDict = {}
2142 for m, mapElement in enumerate(mapElements):
2143 mapType = wcsf.mapCollection.getMapType(mapElement)
2144 if mapType == "Color":
2145 # DCR fit should not go into the generic WCS.
2146 continue
2147 mapDict[mapElement] = {"Type": mapType}
2149 if mapType == "Poly":
2150 mapCoefficients = mapParams[mapElement]
2151 mapDict[mapElement]["Coefficients"] = mapCoefficients
2153 # The RA and Dec of the visit are needed for the last step of
2154 # the mapping from the visit tangent plane to RA and Dec
2155 outWCS = self._make_afw_wcs(
2156 mapDict,
2157 exposureInfo.ras[v] * lsst.geom.radians,
2158 exposureInfo.decs[v] * lsst.geom.radians,
2159 doNormalizePixels=True,
2160 xScale=xscale,
2161 yScale=yscale,
2162 )
2164 catalog[d].setId(detector)
2165 catalog[d].setWcs(outWCS)
2166 catalog.sort()
2167 catalogs[visit] = catalog
2168 if self.config.useColor:
2169 colorVisits = np.array(list(colorFits.keys()))
2170 colorRA = np.array([colorFits[vis][0] for vis in colorVisits])
2171 colorDec = np.array([colorFits[vis][1] for vis in colorVisits])
2172 colorFits = {"visit": colorVisits, "raCoefficient": colorRA, "decCoefficient": colorDec}
2174 return catalogs, cameraParams, colorFits, camera, partialOutputs
2176 def _compute_model_params(self, wcsf):
2177 """Get the WCS model parameters and covariance and convert to a
2178 dictionary that will be readable as a pandas dataframe or other table.
2180 Parameters
2181 ----------
2182 wcsf : `wcsfit.WCSFit`
2183 WCSFit object, assumed to have fit model.
2185 Returns
2186 -------
2187 modelParams : `dict`
2188 Parameters and covariance of the best-fit WCS model.
2189 """
2190 modelParamDict = wcsf.mapCollection.getParamDict()
2191 modelCovariance = wcsf.getModelCovariance()
2193 modelParams = {k: [] for k in ["mapName", "coordinate", "parameter", "coefficientNumber"]}
2194 i = 0
2195 for mapName, params in modelParamDict.items():
2196 nCoeffs = len(params)
2197 # There are an equal number of x and y coordinate parameters
2198 nCoordCoeffs = nCoeffs // 2
2199 modelParams["mapName"].extend([mapName] * nCoeffs)
2200 modelParams["coordinate"].extend(["x"] * nCoordCoeffs)
2201 modelParams["coordinate"].extend(["y"] * nCoordCoeffs)
2202 modelParams["parameter"].extend(params)
2203 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2204 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2206 for p in range(nCoeffs):
2207 if p < nCoordCoeffs:
2208 coord = "x"
2209 else:
2210 coord = "y"
2211 modelParams[f"{mapName}_{coord}_{p}_cov"] = modelCovariance[i]
2212 i += 1
2214 # Convert the dictionary values from lists to numpy arrays.
2215 for key, value in modelParams.items():
2216 modelParams[key] = np.array(value)
2218 return modelParams
2221class GbdesAstrometricMultibandFitConnections(
2222 GbdesAstrometricFitConnections, dimensions=("skymap", "tract", "instrument")
2223):
2224 outputCatalog = pipeBase.connectionTypes.Output(
2225 doc=(
2226 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2227 "plane coordinates and chisq values."
2228 ),
2229 name="gbdesAstrometricMultibandFit_fitStars",
2230 storageClass="ArrowNumpyDict",
2231 dimensions=("instrument", "skymap", "tract"),
2232 )
2233 starCatalog = pipeBase.connectionTypes.Output(
2234 doc=(
2235 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2236 "fitProperMotion is True."
2237 ),
2238 name="gbdesAstrometricMultibandFit_starCatalog",
2239 storageClass="ArrowNumpyDict",
2240 dimensions=("instrument", "skymap", "tract"),
2241 )
2242 modelParams = pipeBase.connectionTypes.Output(
2243 doc="WCS parameters and covariance.",
2244 name="gbdesAstrometricMultibandFit_modelParams",
2245 storageClass="ArrowNumpyDict",
2246 dimensions=("instrument", "skymap", "tract"),
2247 )
2250class GbdesAstrometricMultibandFitConfig(
2251 GbdesAstrometricFitConfig, pipelineConnections=GbdesAstrometricMultibandFitConnections
2252):
2253 pass
2256class GbdesAstrometricMultibandFitTask(GbdesAstrometricFitTask):
2257 """Calibrate the WCS across multiple visits in multiple filters of the same
2258 field using the GBDES package.
2259 """
2261 ConfigClass = GbdesAstrometricMultibandFitConfig
2262 _DefaultName = "gbdesAstrometricMultibandFit"
2265class GbdesGlobalAstrometricFitConnections(
2266 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter")
2267):
2268 inputVisitSummaries = pipeBase.connectionTypes.Input(
2269 doc=(
2270 "Per-visit consolidated exposure metadata built from calexps. "
2271 "These catalogs use detector id for the id and must be sorted for "
2272 "fast lookups of a detector."
2273 ),
2274 name="visitSummary",
2275 storageClass="ExposureCatalog",
2276 dimensions=("instrument", "visit"),
2277 multiple=True,
2278 )
2279 referenceCatalog = pipeBase.connectionTypes.PrerequisiteInput(
2280 doc="The astrometry reference catalog to match to loaded input catalog sources.",
2281 name="gaia_dr3_20230707",
2282 storageClass="SimpleCatalog",
2283 dimensions=("skypix",),
2284 deferLoad=True,
2285 multiple=True,
2286 )
2287 colorCatalog = pipeBase.connectionTypes.Input(
2288 doc="The catalog of magnitudes to match to input sources.",
2289 name="fgcm_Cycle4_StandardStars",
2290 storageClass="SimpleCatalog",
2291 dimensions=("instrument",),
2292 )
2293 isolatedStarSources = pipeBase.connectionTypes.Input(
2294 doc="Catalog of matched sources.",
2295 name="isolated_star_presources",
2296 storageClass="DataFrame",
2297 dimensions=(
2298 "instrument",
2299 "skymap",
2300 "tract",
2301 ),
2302 multiple=True,
2303 deferLoad=True,
2304 )
2305 isolatedStarCatalogs = pipeBase.connectionTypes.Input(
2306 doc="Catalog of objects corresponding to the isolatedStarSources.",
2307 name="isolated_star_presource_associations",
2308 storageClass="DataFrame",
2309 dimensions=(
2310 "instrument",
2311 "skymap",
2312 "tract",
2313 ),
2314 multiple=True,
2315 deferLoad=True,
2316 )
2317 inputCameraModel = pipeBase.connectionTypes.PrerequisiteInput(
2318 doc="Camera parameters to use for 'device' part of model",
2319 name="gbdesAstrometricFit_cameraModel",
2320 storageClass="ArrowNumpyDict",
2321 dimensions=("instrument", "physical_filter"),
2322 )
2323 inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
2324 doc="Input camera to use when constructing camera from astrometric model.",
2325 name="camera",
2326 storageClass="Camera",
2327 dimensions=("instrument",),
2328 isCalibration=True,
2329 )
2330 outputWcs = pipeBase.connectionTypes.Output(
2331 doc=(
2332 "Per-visit world coordinate systems derived from the fitted model. These catalogs only contain "
2333 "entries for detectors with an output, and use the detector id for the catalog id, sorted on id "
2334 "for fast lookups of a detector."
2335 ),
2336 name="gbdesGlobalAstrometricFitSkyWcsCatalog",
2337 storageClass="ExposureCatalog",
2338 dimensions=("instrument", "visit"),
2339 multiple=True,
2340 )
2341 outputCatalog = pipeBase.connectionTypes.Output(
2342 doc=(
2343 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2344 "plane coordinates and chisq values."
2345 ),
2346 name="gbdesGlobalAstrometricFit_fitStars",
2347 storageClass="ArrowNumpyDict",
2348 dimensions=("instrument", "physical_filter"),
2349 )
2350 starCatalog = pipeBase.connectionTypes.Output(
2351 doc=(
2352 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2353 "fitProperMotion is True."
2354 ),
2355 name="gbdesGlobalAstrometricFit_starCatalog",
2356 storageClass="ArrowNumpyDict",
2357 dimensions=("instrument", "physical_filter"),
2358 )
2359 modelParams = pipeBase.connectionTypes.Output(
2360 doc="WCS parameters and covariance.",
2361 name="gbdesGlobalAstrometricFit_modelParams",
2362 storageClass="ArrowNumpyDict",
2363 dimensions=("instrument", "physical_filter"),
2364 )
2365 outputCameraModel = pipeBase.connectionTypes.Output(
2366 doc="Camera parameters to use for 'device' part of model",
2367 name="gbdesAstrometricFit_cameraModel",
2368 storageClass="ArrowNumpyDict",
2369 dimensions=("instrument", "physical_filter"),
2370 )
2371 dcrCoefficients = pipeBase.connectionTypes.Output(
2372 doc="Per-visit coefficients for DCR correction.",
2373 name="gbdesGlobalAstrometricFit_dcrCoefficients",
2374 storageClass="ArrowNumpyDict",
2375 )
2376 camera = pipeBase.connectionTypes.Output(
2377 doc="Camera object constructed using the per-detector part of the astrometric model",
2378 name="gbdesGlobalAstrometricFitCamera",
2379 storageClass="Camera",
2380 dimensions=("instrument", "physical_filter"),
2381 )
2383 def getSpatialBoundsConnections(self):
2384 return ("inputVisitSummaries",)
2386 def __init__(self, *, config=None):
2387 super().__init__(config=config)
2389 if not self.config.useColor:
2390 self.inputs.remove("colorCatalog")
2391 self.outputs.remove("dcrCoefficients")
2392 if not self.config.saveModelParams:
2393 self.outputs.remove("modelParams")
2394 if not self.config.useInputCameraModel:
2395 self.prerequisiteInputs.remove("inputCameraModel")
2396 if not self.config.saveCameraModel:
2397 self.outputs.remove("outputCameraModel")
2398 if not self.config.saveCameraObject:
2399 self.prerequisiteInputs.remove("inputCamera")
2400 self.outputs.remove("camera")
2403class GbdesGlobalAstrometricFitConfig(
2404 GbdesAstrometricFitConfig, pipelineConnections=GbdesGlobalAstrometricFitConnections
2405):
2406 visitOverlap = pexConfig.Field(
2407 dtype=float,
2408 default=1.0,
2409 doc=(
2410 "The linkage distance threshold above which clustered groups of visits will not be merged "
2411 "together in an agglomerative clustering algorithm. The linkage distance is calculated using the "
2412 "minimum distance between the field-of-view centers of a given visit and all other visits in a "
2413 "group, and is in units of the field-of-view radius. The resulting groups of visits define the "
2414 "fields for the astrometric fit."
2415 ),
2416 )
2419class GbdesGlobalAstrometricFitTask(GbdesAstrometricFitTask):
2420 """Calibrate the WCS across multiple visits and multiple fields using the
2421 GBDES package.
2423 This class assumes that the input visits can be separated into contiguous
2424 groups, for which an individual group covers an area of less than a
2425 hemisphere.
2426 """
2428 ConfigClass = GbdesGlobalAstrometricFitConfig
2429 _DefaultName = "gbdesAstrometricFit"
2431 def runQuantum(self, butlerQC, inputRefs, outputRefs):
2432 # We override runQuantum to set up the refObjLoaders
2433 inputs = butlerQC.get(inputRefs)
2435 instrumentName = butlerQC.quantum.dataId["instrument"]
2437 # Ensure the inputs are in a consistent and deterministic order
2438 inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]])
2439 inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()]
2440 inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog])
2441 inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()]
2442 inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]])
2443 inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()]
2444 inputIsolatedStarSourceTracts = np.array(
2445 [isolatedStarSource.dataId["tract"] for isolatedStarSource in inputs["isolatedStarSources"]]
2446 )
2447 inputIsolatedStarCatalogTracts = np.array(
2448 [isolatedStarCatalog.dataId["tract"] for isolatedStarCatalog in inputs["isolatedStarCatalogs"]]
2449 )
2450 for tract in inputIsolatedStarCatalogTracts:
2451 if tract not in inputIsolatedStarSourceTracts:
2452 raise RuntimeError(f"tract {tract} in isolated_star_cats but not isolated_star_sources")
2453 inputs["isolatedStarSources"] = np.array(
2454 [inputs["isolatedStarSources"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2455 )
2456 inputs["isolatedStarCatalogs"] = np.array(
2457 [inputs["isolatedStarCatalogs"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2458 )
2460 refConfig = LoadReferenceObjectsConfig()
2461 if self.config.applyRefCatProperMotion:
2462 refConfig.requireProperMotion = True
2463 refObjectLoader = ReferenceObjectLoader(
2464 dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs],
2465 refCats=inputs.pop("referenceCatalog"),
2466 config=refConfig,
2467 log=self.log,
2468 )
2469 if self.config.useColor:
2470 colorCatalog = inputs.pop("colorCatalog")
2471 else:
2472 colorCatalog = None
2474 try:
2475 output = self.run(
2476 **inputs,
2477 instrumentName=instrumentName,
2478 refObjectLoader=refObjectLoader,
2479 colorCatalog=colorCatalog,
2480 )
2481 except pipeBase.AlgorithmError as e:
2482 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2483 # No partial outputs for butler to put
2484 raise error from e
2486 for outputRef in outputRefs.outputWcs:
2487 visit = outputRef.dataId["visit"]
2488 butlerQC.put(output.outputWcss[visit], outputRef)
2489 butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
2490 butlerQC.put(output.starCatalog, outputRefs.starCatalog)
2491 if self.config.saveModelParams:
2492 butlerQC.put(output.modelParams, outputRefs.modelParams)
2493 if self.config.saveCameraModel:
2494 butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel)
2495 if self.config.saveCameraObject:
2496 butlerQC.put(output.camera, outputRefs.camera)
2497 if self.config.useColor:
2498 butlerQC.put(output.colorParams, outputRefs.dcrCoefficients)
2499 if output.partialOutputs:
2500 e = RuntimeError("Some visits were dropped because data was insufficient to fit model.")
2501 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2502 raise error from e
2504 def run(
2505 self,
2506 inputVisitSummaries,
2507 isolatedStarSources,
2508 isolatedStarCatalogs,
2509 instrumentName="",
2510 refEpoch=None,
2511 refObjectLoader=None,
2512 inputCameraModel=None,
2513 colorCatalog=None,
2514 inputCamera=None,
2515 ):
2516 """Run the WCS fit for a given set of visits
2518 Parameters
2519 ----------
2520 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
2521 List of catalogs with per-detector summary information.
2522 isolatedStarSources : `list` [`DeferredDatasetHandle`]
2523 List of handles pointing to isolated star sources.
2524 isolatedStarCatalog: `list` [`DeferredDatasetHandle`]
2525 List of handles pointing to isolated star catalogs.
2526 instrumentName : `str`, optional
2527 Name of the instrument used. This is only used for labelling.
2528 refEpoch : `float`, optional
2529 Epoch of the reference objects in MJD.
2530 refObjectLoader : instance of
2531 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`,
2532 optional
2533 Reference object loader instance.
2534 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
2535 Parameters to use for the device part of the model.
2536 colorCatalog : `lsst.afw.table.SimpleCatalog`
2537 Catalog containing object coordinates and magnitudes.
2538 inputCamera : `lsst.afw.cameraGeom.Camera`, optional
2539 Camera to be used as template when constructing new camera.
2541 Returns
2542 -------
2543 result : `lsst.pipe.base.Struct`
2544 ``outputWcss`` : `list` [`lsst.afw.table.ExposureCatalog`]
2545 List of exposure catalogs (one per visit) with the WCS for each
2546 detector set by the new fitted WCS.
2547 ``fitModel`` : `wcsfit.WCSFit`
2548 Model-fitting object with final model parameters.
2549 ``outputCatalog`` : `pyarrow.Table`
2550 Catalog with fit residuals of all sources used.
2551 ``starCatalog`` : `pyarrow.Table`
2552 Catalog with best-fit positions of the objects fit.
2553 ``modelParams`` : `dict`
2554 Parameters and covariance of the best-fit WCS model.
2555 ``cameraModelParams`` : `dict` [`str`, `np.ndarray`]
2556 Parameters of the device part of the model, in the format
2557 needed as input for future runs.
2558 ``colorParams`` : `dict` [`int`, `np.ndarray`]
2559 DCR parameters fit in RA and Dec directions for each visit.
2560 ``camera`` : `lsst.afw.cameraGeom.Camera`
2561 Camera object constructed from the per-detector model.
2562 """
2563 self.log.info("Gather instrument, exposure, and field info")
2565 # Get information about the extent of the input visits
2566 fields, fieldRegions = self._prep_sky(inputVisitSummaries)
2568 # Get RA, Dec, MJD, etc., for the input visits
2569 exposureInfo, exposuresHelper, extensionInfo, instruments = self._get_exposure_info(
2570 inputVisitSummaries, fieldRegions=fieldRegions
2571 )
2573 self.log.info("Load associated sources")
2574 medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="jyear").mjd
2575 allRefObjects, allRefCovariances = {}, {}
2576 for f, fieldRegion in fieldRegions.items():
2577 refObjects, refCovariance = self._load_refcat(
2578 refObjectLoader, extensionInfo, epoch=medianEpoch, region=fieldRegion
2579 )
2580 allRefObjects[f] = refObjects
2581 allRefCovariances[f] = refCovariance
2583 associations, sourceDict = self._associate_from_isolated_sources(
2584 isolatedStarSources, isolatedStarCatalogs, extensionInfo, allRefObjects
2585 )
2587 self.log.info("Fit the WCSs")
2588 # Set up a YAML-type string using the config variables and a sample
2589 # visit
2590 inputYaml, mapTemplate = self.make_yaml(
2591 inputVisitSummaries[0],
2592 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
2593 )
2595 # Set the verbosity level for WCSFit from the task log level.
2596 # TODO: DM-36850, Add lsst.log to gbdes so that log messages are
2597 # properly propagated.
2598 loglevel = self.log.getEffectiveLevel()
2599 if loglevel >= self.log.WARNING:
2600 verbose = 0
2601 elif loglevel == self.log.INFO:
2602 verbose = 1
2603 else:
2604 verbose = 2
2606 # Set up the WCS-fitting class using the source matches from the
2607 # isolated star sources plus the reference catalog.
2608 wcsf = wcsfit.WCSFit(
2609 fields,
2610 instruments,
2611 exposuresHelper,
2612 extensionInfo.visitIndex,
2613 extensionInfo.detectorIndex,
2614 inputYaml,
2615 extensionInfo.wcs,
2616 associations.sequence,
2617 associations.extn,
2618 associations.obj,
2619 sysErr=self.config.systematicError,
2620 refSysErr=self.config.referenceSystematicError,
2621 usePM=self.config.fitProperMotion,
2622 verbose=verbose,
2623 )
2625 # Add the science and reference sources
2626 self._add_objects(wcsf, sourceDict, extensionInfo)
2627 for f in fieldRegions.keys():
2628 self._add_ref_objects(
2629 wcsf, allRefObjects[f], allRefCovariances[f], extensionInfo, fieldIndex=-1 * f
2630 )
2631 if self.config.useColor:
2632 self._add_color_objects(wcsf, colorCatalog)
2634 # Do the WCS fit
2635 try:
2636 wcsf.fit(
2637 reserveFraction=self.config.fitReserveFraction,
2638 randomNumberSeed=self.config.fitReserveRandomSeed,
2639 clipThresh=self.config.clipThresh,
2640 clipFraction=self.config.clipFraction,
2641 )
2642 except RuntimeError as e:
2643 if "Cholesky decomposition failed" in str(e):
2644 raise CholeskyError() from e
2645 else:
2646 raise
2648 self.log.info("WCS fitting done")
2650 outputWcss, cameraParams, colorParams, camera, partialOutputs = self._make_outputs(
2651 wcsf,
2652 inputVisitSummaries,
2653 exposureInfo,
2654 mapTemplate,
2655 extensionInfo,
2656 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
2657 inputCamera=(inputCamera if self.config.buildCamera else None),
2658 )
2659 outputCatalog = wcsf.getOutputCatalog()
2660 outputCatalog["exposureName"] = np.array(outputCatalog["exposureName"])
2661 outputCatalog["deviceName"] = np.array(outputCatalog["deviceName"])
2663 starCatalog = wcsf.getStarCatalog()
2664 modelParams = self._compute_model_params(wcsf) if self.config.saveModelParams else None
2666 return pipeBase.Struct(
2667 outputWcss=outputWcss,
2668 fitModel=wcsf,
2669 outputCatalog=outputCatalog,
2670 starCatalog=starCatalog,
2671 modelParams=modelParams,
2672 cameraModelParams=cameraParams,
2673 colorParams=colorParams,
2674 camera=camera,
2675 partialOutputs=partialOutputs,
2676 )
2678 def _prep_sky(self, inputVisitSummaries):
2679 """Cluster the input visits into disjoint groups that will define
2680 separate fields in the astrometric fit, and, for each group, get the
2681 convex hull around all of its component visits.
2683 The groups are created such that each visit overlaps with at least one
2684 other visit in the same group by the `visitOverlap` amount, which is
2685 calculated as a fraction of the field-of-view radius, and no visits
2686 from separate groups overlap by more than this amount.
2688 Paramaters
2689 ----------
2690 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
2691 List of catalogs with per-detector summary information.
2693 Returns
2694 -------
2695 fields : `wcsfit.Fields`
2696 Object with field information.
2697 fieldRegions : `dict` [`int`, `lsst.sphgeom.ConvexPolygon`]
2698 Dictionary of regions encompassing each group of input visits,
2699 keyed by an arbitrary index.
2700 """
2701 allDetectorCorners = []
2702 mjds = []
2703 radecs = []
2704 radii = []
2705 for visSum in inputVisitSummaries:
2706 detectorCorners = [
2707 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
2708 for (ra, dec) in zip(visSum["raCorners"].ravel(), visSum["decCorners"].ravel())
2709 if (np.isfinite(ra) and (np.isfinite(dec)))
2710 ]
2711 if len(detectorCorners) == 0:
2712 # Skip this visit if none of the detectors have finite ra/dec
2713 # corners, which happens when the WCSs are missing. The visit
2714 # will get formally dropped elsewhere.
2715 continue
2716 allDetectorCorners.append(detectorCorners)
2718 # Get center and approximate radius of field of view
2719 boundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(detectorCorners).getBoundingCircle()
2720 center = lsst.geom.SpherePoint(boundingCircle.getCenter())
2721 ra = center.getRa().asDegrees()
2722 dec = center.getDec().asDegrees()
2723 radecs.append([ra, dec])
2724 radius = boundingCircle.getOpeningAngle()
2725 radii.append(radius)
2727 obsDate = visSum[0].getVisitInfo().getDate()
2728 obsMJD = obsDate.get(obsDate.MJD)
2729 mjds.append(obsMJD)
2731 # Find groups of visits where any one of the visits overlaps another by
2732 # a given fraction of the field-of-view radius.
2733 distance = self.config.visitOverlap * np.median(radii)
2734 clustering = AgglomerativeClustering(
2735 distance_threshold=distance.asDegrees(), n_clusters=None, linkage="single"
2736 )
2737 clusters = clustering.fit(np.array(radecs))
2739 medianMJD = np.median(mjds)
2740 medianEpoch = astropy.time.Time(medianMJD, format="mjd").jyear
2742 fieldNames = []
2743 fieldRAs = []
2744 fieldDecs = []
2745 epochs = []
2746 fieldRegions = {}
2748 for i in range(clusters.n_clusters_):
2749 fieldInd = clusters.labels_ == i
2750 # Concatenate the lists of all detector corners that are in this
2751 # field
2752 fieldDetectors = sum([allDetectorCorners[f] for f, fInd in enumerate(fieldInd) if fInd], [])
2753 hull = lsst.sphgeom.ConvexPolygon.convexHull(fieldDetectors)
2754 center = lsst.geom.SpherePoint(hull.getCentroid())
2755 ra = center.getRa().asDegrees()
2756 dec = center.getDec().asDegrees()
2758 fieldRegions[i] = hull
2759 fieldNames.append(str(i))
2760 fieldRAs.append(ra)
2761 fieldDecs.append(dec)
2762 # Use the same median epoch for all fields so that the final object
2763 # positions are calculated for the same epoch.
2764 epochs.append(medianEpoch)
2766 fields = wcsfit.Fields(fieldNames, fieldRAs, fieldDecs, epochs)
2768 return fields, fieldRegions
2770 def _associate_from_isolated_sources(
2771 self, isolatedStarSourceRefs, isolatedStarCatalogRefs, extensionInfo, refObjects
2772 ):
2773 """Match the input catalog of isolated stars with the reference catalog
2774 and transform the combined isolated star sources and reference source
2775 into the format needed for gbdes.
2777 Parameters
2778 ----------
2779 isolatedStarSourceRefs : `list` [`DeferredDatasetHandle`]
2780 List of handles pointing to isolated star sources.
2781 isolatedStarCatalogRefs: `list` [`DeferredDatasetHandle`]
2782 List of handles pointing to isolated star catalogs.
2783 extensionInfo : `lsst.pipe.base.Struct`
2784 Struct containing properties for each extension (visit/detector).
2785 ``visit`` : `np.ndarray`
2786 Name of the visit for this extension.
2787 ``detector`` : `np.ndarray`
2788 Name of the detector for this extension.
2789 ``visitIndex` : `np.ndarray` [`int`]
2790 Index of visit for this extension.
2791 ``detectorIndex`` : `np.ndarray` [`int`]
2792 Index of the detector for this extension.
2793 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2794 Initial WCS for this extension.
2795 ``extensionType`` : `np.ndarray` [`str`]
2796 "SCIENCE" or "REFERENCE".
2797 refObjects : `dict`
2798 Dictionary of dictionaries containing the position and error
2799 information of reference objects.
2801 Returns
2802 -------
2803 associations : `lsst.pipe.base.Struct`
2804 Struct containing the associations of sources with objects.
2805 sourceDict : `dict` [`int`, [`int`, [`str`, `list` [`float`]]]]
2806 Dictionary containing the source centroids for each visit.
2807 """
2808 sequences = []
2809 extensions = []
2810 object_indices = []
2812 sourceColumns = ["x", "y", "xErr", "yErr", "ixx", "ixy", "iyy", "obj_index", "visit", "detector"]
2813 catalogColumns = ["ra", "dec"]
2815 sourceDict = dict([(visit, {}) for visit in np.unique(extensionInfo.visit)])
2816 for visit, detector in zip(extensionInfo.visit, extensionInfo.detector):
2817 sourceDict[visit][detector] = {"x": [], "y": [], "xCov": [], "yCov": [], "xyCov": []}
2819 for isolatedStarCatalogRef, isolatedStarSourceRef in zip(
2820 isolatedStarCatalogRefs, isolatedStarSourceRefs
2821 ):
2822 isolatedStarCatalog = isolatedStarCatalogRef.get(parameters={"columns": catalogColumns})
2823 isolatedStarSources = isolatedStarSourceRef.get(parameters={"columns": sourceColumns})
2824 if len(isolatedStarCatalog) == 0:
2825 # This is expected when only one visit overlaps with a given
2826 # tract, meaning that no sources can be associated.
2827 self.log.debug(
2828 "Skipping tract %d, which has no associated isolated stars",
2829 isolatedStarCatalogRef.dataId["tract"],
2830 )
2831 continue
2833 # Match the reference stars to the existing isolated stars, then
2834 # insert the reference stars into the isolated star sources.
2835 allVisits = np.copy(isolatedStarSources["visit"])
2836 allDetectors = np.copy(isolatedStarSources["detector"])
2837 allObjectIndices = np.copy(isolatedStarSources["obj_index"])
2838 issIndices = np.copy(isolatedStarSources.index)
2839 for f, regionRefObjects in refObjects.items():
2840 # Use the same matching technique that is done in
2841 # isolatedStarAssociation and fgcmBuildFromIsolatedStars.
2842 with Matcher(
2843 isolatedStarCatalog["ra"].to_numpy(), isolatedStarCatalog["dec"].to_numpy()
2844 ) as matcher:
2845 idx, idx_isoStarCat, idx_refObjects, d = matcher.query_radius(
2846 np.array(regionRefObjects["ra"]),
2847 np.array(regionRefObjects["dec"]),
2848 self.config.matchRadius / 3600.0,
2849 return_indices=True,
2850 )
2852 refSort = np.searchsorted(isolatedStarSources["obj_index"], idx_isoStarCat)
2853 refDetector = np.ones(len(idx_isoStarCat)) * -1
2854 # The "visit" for the reference catalogs is the field times -1.
2855 refVisit = np.ones(len(idx_isoStarCat)) * f * -1
2857 allVisits = np.insert(allVisits, refSort, refVisit)
2858 allDetectors = np.insert(allDetectors, refSort, refDetector)
2859 allObjectIndices = np.insert(allObjectIndices, refSort, idx_isoStarCat)
2860 issIndices = np.insert(issIndices, refSort, idx_refObjects)
2862 # Loop through the associated sources to convert them to the gbdes
2863 # format, which requires the extension index, the source's index in
2864 # the input table, and a sequence number corresponding to the
2865 # object with which it is associated.
2866 sequence = 0
2867 obj_index = allObjectIndices[0]
2868 for visit, detector, row, obj_ind in zip(allVisits, allDetectors, issIndices, allObjectIndices):
2869 extensionIndex = np.flatnonzero(
2870 (extensionInfo.visit == visit) & (extensionInfo.detector == detector)
2871 )
2872 if len(extensionIndex) == 0:
2873 # This happens for runs where you are not using all the
2874 # visits overlapping a given tract that were included in
2875 # the isolated star association task."
2876 continue
2877 else:
2878 extensionIndex = extensionIndex[0]
2880 extensions.append(extensionIndex)
2881 if visit <= 0:
2882 object_indices.append(row)
2883 else:
2884 object_indices.append(len(sourceDict[visit][detector]["x"]))
2885 source = isolatedStarSources.loc[row]
2886 sourceDict[visit][detector]["x"].append(source["x"])
2887 sourceDict[visit][detector]["y"].append(source["y"])
2888 xCov = source["xErr"] ** 2
2889 yCov = source["yErr"] ** 2
2890 xyCov = source["ixy"] * (xCov + yCov) / (source["ixx"] + source["iyy"])
2891 # TODO: add correct xyErr if DM-7101 is ever done.
2892 sourceDict[visit][detector]["xCov"].append(xCov)
2893 sourceDict[visit][detector]["yCov"].append(yCov)
2894 sourceDict[visit][detector]["xyCov"].append(xyCov)
2895 if obj_ind != obj_index:
2896 sequence = 0
2897 sequences.append(sequence)
2898 obj_index = obj_ind
2899 sequence += 1
2900 else:
2901 sequences.append(sequence)
2902 sequence += 1
2904 associations = pipeBase.Struct(extn=extensions, obj=object_indices, sequence=sequences)
2905 return associations, sourceDict
2907 def _add_objects(self, wcsf, sourceDict, extensionInfo):
2908 """Add science sources to the wcsfit.WCSFit object.
2910 Parameters
2911 ----------
2912 wcsf : `wcsfit.WCSFit`
2913 WCS-fitting object.
2914 sourceDict : `dict`
2915 Dictionary containing the source centroids for each visit.
2916 extensionInfo : `lsst.pipe.base.Struct`
2917 Struct containing properties for each extension (visit/detector).
2918 ``visit`` : `np.ndarray`
2919 Name of the visit for this extension.
2920 ``detector`` : `np.ndarray`
2921 Name of the detector for this extension.
2922 ``visitIndex` : `np.ndarray` [`int`]
2923 Index of visit for this extension.
2924 ``detectorIndex`` : `np.ndarray` [`int`]
2925 Index of the detector for this extension.
2926 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2927 Initial WCS for this extension.
2928 ``extensionType`` : `np.ndarray` [`str`]
2929 "SCIENCE" or "REFERENCE".
2930 """
2931 for visit, visitSources in sourceDict.items():
2932 # Visit numbers equal or below zero connote the reference catalog.
2933 if visit <= 0:
2934 # This "visit" number corresponds to a reference catalog.
2935 continue
2937 for detector, sourceCat in visitSources.items():
2938 extensionIndex = np.flatnonzero(
2939 (extensionInfo.visit == visit) & (extensionInfo.detector == detector)
2940 )[0]
2942 d = {
2943 "x": np.array(sourceCat["x"]),
2944 "y": np.array(sourceCat["y"]),
2945 "xCov": np.array(sourceCat["xCov"]),
2946 "yCov": np.array(sourceCat["yCov"]),
2947 "xyCov": np.array(sourceCat["xyCov"]),
2948 }
2949 wcsf.setObjects(
2950 extensionIndex,
2951 d,
2952 "x",
2953 "y",
2954 ["xCov", "yCov", "xyCov"],
2955 defaultColor=self.config.referenceColor,
2956 )
2959class GbdesGlobalAstrometricMultibandFitConnections(
2960 GbdesGlobalAstrometricFitConnections,
2961 dimensions=("instrument",),
2962):
2963 outputCatalog = pipeBase.connectionTypes.Output(
2964 doc=(
2965 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2966 "plane coordinates and chisq values."
2967 ),
2968 name="gbdesGlobalAstrometricMultibandFit_fitStars",
2969 storageClass="ArrowNumpyDict",
2970 dimensions=("instrument",),
2971 )
2972 starCatalog = pipeBase.connectionTypes.Output(
2973 doc=(
2974 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2975 "fitProperMotion is True."
2976 ),
2977 name="gbdesGlobalAstrometricMultibandFit_starCatalog",
2978 storageClass="ArrowNumpyDict",
2979 dimensions=("instrument",),
2980 )
2981 modelParams = pipeBase.connectionTypes.Output(
2982 doc="WCS parameters and covariance.",
2983 name="gbdesGlobalAstrometricMultibandFit_modelParams",
2984 storageClass="ArrowNumpyDict",
2985 dimensions=("instrument",),
2986 )
2989class GbdesGlobalAstrometricMultibandFitConfig(
2990 GbdesAstrometricFitConfig,
2991 pipelineConnections=GbdesGlobalAstrometricMultibandFitConnections,
2992):
2993 """Configuration for the GbdesGlobalAstrometricMultibandFitTask"""
2995 pass
2998class GbdesGlobalAstrometricMultibandFitTask(GbdesGlobalAstrometricFitTask):
2999 """Calibrate the WCS across multiple visits in multiple filters and
3000 multiple fields using the GBDES package.
3001 """
3003 ConfigClass = GbdesGlobalAstrometricMultibandFitConfig
3004 _DefaultName = "gbdesAstrometricMultibandFit"