Coverage for python / lsst / drp / tasks / gbdesAstrometricFit.py: 11%
1027 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-15 00:19 +0000
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 # Make dictionary of bboxes for each detector.
2089 detectorBBoxes = {}
2090 for detector in exposureInfo.detectors:
2091 for visitSummary in visitSummaryTables:
2092 if detInfo := visitSummary.find(detector):
2093 detectorBBoxes[detector] = detInfo.getBBox()
2094 break
2096 catalogs = {}
2097 colorFits = {}
2098 partialOutputs = False
2099 for v, visitSummary in enumerate(visitSummaryTables):
2100 visit = visitSummary[0]["visit"]
2101 if visit not in extensionInfo.visit:
2102 self.log.warning("Visit %d was dropped because no detectors had valid WCSs.", visit)
2103 partialOutputs = True
2104 continue
2106 visitMaps = wcsf.mapCollection.orderAtoms(f"{visit}")
2107 if self.config.useColor:
2108 colorMap = visitMaps.pop(visitMaps.index(f"{visit}/dcr"))
2109 colorFits[visit] = mapParams[colorMap]
2110 visitMap = visitMaps[0]
2111 visitMapType = wcsf.mapCollection.getMapType(visitMap)
2112 if (visitMap not in mapParams) and (visitMapType != "Identity"):
2113 self.log.warning("Visit %d was dropped because of an insufficient amount of data.", visit)
2114 partialOutputs = True
2115 continue
2117 catalog = lsst.afw.table.ExposureCatalog(schema)
2118 catalog.resize(len(exposureInfo.detectors))
2119 catalog["visit"] = visit
2121 for d, detector in enumerate(exposureInfo.detectors):
2122 mapName = f"{visit}/{detector}"
2123 if mapName in wcsf.mapCollection.allMapNames():
2124 mapElements = wcsf.mapCollection.orderAtoms(f"{mapName}/base")
2125 catalog[d]["recoveredWcs"] = False
2126 else:
2127 # This extension was not fit, but the WCS can be recovered
2128 # using the maps fit from sources on other visits but the
2129 # same detector and from sources on other detectors from
2130 # this visit.
2131 genericElements = mapTemplate["EXPOSURE/DEVICE/base"]["Elements"]
2132 mapElements = []
2133 band = visitSummary[0]["physical_filter"]
2135 # Go through the generic map components to build the names
2136 # of the specific maps for this extension.
2137 for component in genericElements:
2138 elements = mapTemplate[component]["Elements"]
2139 for element in elements:
2140 # TODO: DM-42519, gbdes sets the "BAND" to the
2141 # instrument name currently. This will need to be
2142 # disambiguated if we run on multiple bands at
2143 # once.
2144 element = element.replace("BAND", str(band))
2145 element = element.replace("EXPOSURE", str(visit))
2146 element = element.replace("DEVICE", str(detector))
2147 mapElements.append(element)
2148 catalog[d]["recoveredWcs"] = True
2149 mapDict = {}
2150 for m, mapElement in enumerate(mapElements):
2151 mapType = wcsf.mapCollection.getMapType(mapElement)
2152 if mapType == "Color":
2153 # DCR fit should not go into the generic WCS.
2154 continue
2155 mapDict[mapElement] = {"Type": mapType}
2157 if mapType == "Poly":
2158 mapCoefficients = mapParams[mapElement]
2159 mapDict[mapElement]["Coefficients"] = mapCoefficients
2161 # The RA and Dec of the visit are needed for the last step of
2162 # the mapping from the visit tangent plane to RA and Dec
2163 outWCS = self._make_afw_wcs(
2164 mapDict,
2165 exposureInfo.ras[v] * lsst.geom.radians,
2166 exposureInfo.decs[v] * lsst.geom.radians,
2167 doNormalizePixels=True,
2168 xScale=xscale,
2169 yScale=yscale,
2170 )
2172 catalog[d].setId(detector)
2173 catalog[d].setBBox(detectorBBoxes[detector])
2174 catalog[d].setWcs(outWCS)
2175 catalog.sort()
2176 catalogs[visit] = catalog
2177 if self.config.useColor:
2178 colorVisits = np.array(list(colorFits.keys()))
2179 colorRA = np.array([colorFits[vis][0] for vis in colorVisits])
2180 colorDec = np.array([colorFits[vis][1] for vis in colorVisits])
2181 colorFits = {"visit": colorVisits, "raCoefficient": colorRA, "decCoefficient": colorDec}
2183 return catalogs, cameraParams, colorFits, camera, partialOutputs
2185 def _compute_model_params(self, wcsf):
2186 """Get the WCS model parameters and covariance and convert to a
2187 dictionary that will be readable as a pandas dataframe or other table.
2189 Parameters
2190 ----------
2191 wcsf : `wcsfit.WCSFit`
2192 WCSFit object, assumed to have fit model.
2194 Returns
2195 -------
2196 modelParams : `dict`
2197 Parameters and covariance of the best-fit WCS model.
2198 """
2199 modelParamDict = wcsf.mapCollection.getParamDict()
2200 modelCovariance = wcsf.getModelCovariance()
2202 modelParams = {k: [] for k in ["mapName", "coordinate", "parameter", "coefficientNumber"]}
2203 i = 0
2204 for mapName, params in modelParamDict.items():
2205 nCoeffs = len(params)
2206 # There are an equal number of x and y coordinate parameters
2207 nCoordCoeffs = nCoeffs // 2
2208 modelParams["mapName"].extend([mapName] * nCoeffs)
2209 modelParams["coordinate"].extend(["x"] * nCoordCoeffs)
2210 modelParams["coordinate"].extend(["y"] * nCoordCoeffs)
2211 modelParams["parameter"].extend(params)
2212 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2213 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2215 for p in range(nCoeffs):
2216 if p < nCoordCoeffs:
2217 coord = "x"
2218 else:
2219 coord = "y"
2220 modelParams[f"{mapName}_{coord}_{p}_cov"] = modelCovariance[i]
2221 i += 1
2223 # Convert the dictionary values from lists to numpy arrays.
2224 for key, value in modelParams.items():
2225 modelParams[key] = np.array(value)
2227 return modelParams
2230class GbdesAstrometricMultibandFitConnections(
2231 GbdesAstrometricFitConnections, dimensions=("skymap", "tract", "instrument")
2232):
2233 outputCatalog = pipeBase.connectionTypes.Output(
2234 doc=(
2235 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2236 "plane coordinates and chisq values."
2237 ),
2238 name="gbdesAstrometricMultibandFit_fitStars",
2239 storageClass="ArrowNumpyDict",
2240 dimensions=("instrument", "skymap", "tract"),
2241 )
2242 starCatalog = pipeBase.connectionTypes.Output(
2243 doc=(
2244 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2245 "fitProperMotion is True."
2246 ),
2247 name="gbdesAstrometricMultibandFit_starCatalog",
2248 storageClass="ArrowNumpyDict",
2249 dimensions=("instrument", "skymap", "tract"),
2250 )
2251 modelParams = pipeBase.connectionTypes.Output(
2252 doc="WCS parameters and covariance.",
2253 name="gbdesAstrometricMultibandFit_modelParams",
2254 storageClass="ArrowNumpyDict",
2255 dimensions=("instrument", "skymap", "tract"),
2256 )
2259class GbdesAstrometricMultibandFitConfig(
2260 GbdesAstrometricFitConfig, pipelineConnections=GbdesAstrometricMultibandFitConnections
2261):
2262 pass
2265class GbdesAstrometricMultibandFitTask(GbdesAstrometricFitTask):
2266 """Calibrate the WCS across multiple visits in multiple filters of the same
2267 field using the GBDES package.
2268 """
2270 ConfigClass = GbdesAstrometricMultibandFitConfig
2271 _DefaultName = "gbdesAstrometricMultibandFit"
2274class GbdesGlobalAstrometricFitConnections(
2275 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter")
2276):
2277 inputVisitSummaries = pipeBase.connectionTypes.Input(
2278 doc=(
2279 "Per-visit consolidated exposure metadata built from calexps. "
2280 "These catalogs use detector id for the id and must be sorted for "
2281 "fast lookups of a detector."
2282 ),
2283 name="visitSummary",
2284 storageClass="ExposureCatalog",
2285 dimensions=("instrument", "visit"),
2286 multiple=True,
2287 )
2288 referenceCatalog = pipeBase.connectionTypes.PrerequisiteInput(
2289 doc="The astrometry reference catalog to match to loaded input catalog sources.",
2290 name="gaia_dr3_20230707",
2291 storageClass="SimpleCatalog",
2292 dimensions=("skypix",),
2293 deferLoad=True,
2294 multiple=True,
2295 )
2296 colorCatalog = pipeBase.connectionTypes.Input(
2297 doc="The catalog of magnitudes to match to input sources.",
2298 name="fgcm_Cycle4_StandardStars",
2299 storageClass="SimpleCatalog",
2300 dimensions=("instrument",),
2301 )
2302 isolatedStarSources = pipeBase.connectionTypes.Input(
2303 doc="Catalog of matched sources.",
2304 name="isolated_star_presources",
2305 storageClass="DataFrame",
2306 dimensions=(
2307 "instrument",
2308 "skymap",
2309 "tract",
2310 ),
2311 multiple=True,
2312 deferLoad=True,
2313 )
2314 isolatedStarCatalogs = pipeBase.connectionTypes.Input(
2315 doc="Catalog of objects corresponding to the isolatedStarSources.",
2316 name="isolated_star_presource_associations",
2317 storageClass="DataFrame",
2318 dimensions=(
2319 "instrument",
2320 "skymap",
2321 "tract",
2322 ),
2323 multiple=True,
2324 deferLoad=True,
2325 )
2326 inputCameraModel = pipeBase.connectionTypes.PrerequisiteInput(
2327 doc="Camera parameters to use for 'device' part of model",
2328 name="gbdesAstrometricFit_cameraModel",
2329 storageClass="ArrowNumpyDict",
2330 dimensions=("instrument", "physical_filter"),
2331 )
2332 inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
2333 doc="Input camera to use when constructing camera from astrometric model.",
2334 name="camera",
2335 storageClass="Camera",
2336 dimensions=("instrument",),
2337 isCalibration=True,
2338 )
2339 outputWcs = pipeBase.connectionTypes.Output(
2340 doc=(
2341 "Per-visit world coordinate systems derived from the fitted model. These catalogs only contain "
2342 "entries for detectors with an output, and use the detector id for the catalog id, sorted on id "
2343 "for fast lookups of a detector."
2344 ),
2345 name="gbdesGlobalAstrometricFitSkyWcsCatalog",
2346 storageClass="ExposureCatalog",
2347 dimensions=("instrument", "visit"),
2348 multiple=True,
2349 )
2350 outputCatalog = pipeBase.connectionTypes.Output(
2351 doc=(
2352 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2353 "plane coordinates and chisq values."
2354 ),
2355 name="gbdesGlobalAstrometricFit_fitStars",
2356 storageClass="ArrowNumpyDict",
2357 dimensions=("instrument", "physical_filter"),
2358 )
2359 starCatalog = pipeBase.connectionTypes.Output(
2360 doc=(
2361 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2362 "fitProperMotion is True."
2363 ),
2364 name="gbdesGlobalAstrometricFit_starCatalog",
2365 storageClass="ArrowNumpyDict",
2366 dimensions=("instrument", "physical_filter"),
2367 )
2368 modelParams = pipeBase.connectionTypes.Output(
2369 doc="WCS parameters and covariance.",
2370 name="gbdesGlobalAstrometricFit_modelParams",
2371 storageClass="ArrowNumpyDict",
2372 dimensions=("instrument", "physical_filter"),
2373 )
2374 outputCameraModel = pipeBase.connectionTypes.Output(
2375 doc="Camera parameters to use for 'device' part of model",
2376 name="gbdesAstrometricFit_cameraModel",
2377 storageClass="ArrowNumpyDict",
2378 dimensions=("instrument", "physical_filter"),
2379 )
2380 dcrCoefficients = pipeBase.connectionTypes.Output(
2381 doc="Per-visit coefficients for DCR correction.",
2382 name="gbdesGlobalAstrometricFit_dcrCoefficients",
2383 storageClass="ArrowNumpyDict",
2384 )
2385 camera = pipeBase.connectionTypes.Output(
2386 doc="Camera object constructed using the per-detector part of the astrometric model",
2387 name="gbdesGlobalAstrometricFitCamera",
2388 storageClass="Camera",
2389 dimensions=("instrument", "physical_filter"),
2390 )
2392 def getSpatialBoundsConnections(self):
2393 return ("inputVisitSummaries",)
2395 def __init__(self, *, config=None):
2396 super().__init__(config=config)
2398 if not self.config.useColor:
2399 self.inputs.remove("colorCatalog")
2400 self.outputs.remove("dcrCoefficients")
2401 if not self.config.saveModelParams:
2402 self.outputs.remove("modelParams")
2403 if not self.config.useInputCameraModel:
2404 self.prerequisiteInputs.remove("inputCameraModel")
2405 if not self.config.saveCameraModel:
2406 self.outputs.remove("outputCameraModel")
2407 if not self.config.saveCameraObject:
2408 self.prerequisiteInputs.remove("inputCamera")
2409 self.outputs.remove("camera")
2412class GbdesGlobalAstrometricFitConfig(
2413 GbdesAstrometricFitConfig, pipelineConnections=GbdesGlobalAstrometricFitConnections
2414):
2415 visitOverlap = pexConfig.Field(
2416 dtype=float,
2417 default=1.0,
2418 doc=(
2419 "The linkage distance threshold above which clustered groups of visits will not be merged "
2420 "together in an agglomerative clustering algorithm. The linkage distance is calculated using the "
2421 "minimum distance between the field-of-view centers of a given visit and all other visits in a "
2422 "group, and is in units of the field-of-view radius. The resulting groups of visits define the "
2423 "fields for the astrometric fit."
2424 ),
2425 )
2428class GbdesGlobalAstrometricFitTask(GbdesAstrometricFitTask):
2429 """Calibrate the WCS across multiple visits and multiple fields using the
2430 GBDES package.
2432 This class assumes that the input visits can be separated into contiguous
2433 groups, for which an individual group covers an area of less than a
2434 hemisphere.
2435 """
2437 ConfigClass = GbdesGlobalAstrometricFitConfig
2438 _DefaultName = "gbdesAstrometricFit"
2440 def runQuantum(self, butlerQC, inputRefs, outputRefs):
2441 # We override runQuantum to set up the refObjLoaders
2442 inputs = butlerQC.get(inputRefs)
2444 instrumentName = butlerQC.quantum.dataId["instrument"]
2446 # Ensure the inputs are in a consistent and deterministic order
2447 inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]])
2448 inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()]
2449 inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog])
2450 inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()]
2451 inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]])
2452 inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()]
2453 inputIsolatedStarSourceTracts = np.array(
2454 [isolatedStarSource.dataId["tract"] for isolatedStarSource in inputs["isolatedStarSources"]]
2455 )
2456 inputIsolatedStarCatalogTracts = np.array(
2457 [isolatedStarCatalog.dataId["tract"] for isolatedStarCatalog in inputs["isolatedStarCatalogs"]]
2458 )
2459 for tract in inputIsolatedStarCatalogTracts:
2460 if tract not in inputIsolatedStarSourceTracts:
2461 raise RuntimeError(f"tract {tract} in isolated_star_cats but not isolated_star_sources")
2462 inputs["isolatedStarSources"] = np.array(
2463 [inputs["isolatedStarSources"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2464 )
2465 inputs["isolatedStarCatalogs"] = np.array(
2466 [inputs["isolatedStarCatalogs"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2467 )
2469 refConfig = LoadReferenceObjectsConfig()
2470 if self.config.applyRefCatProperMotion:
2471 refConfig.requireProperMotion = True
2472 refObjectLoader = ReferenceObjectLoader(
2473 dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs],
2474 refCats=inputs.pop("referenceCatalog"),
2475 config=refConfig,
2476 log=self.log,
2477 )
2478 if self.config.useColor:
2479 colorCatalog = inputs.pop("colorCatalog")
2480 else:
2481 colorCatalog = None
2483 try:
2484 output = self.run(
2485 **inputs,
2486 instrumentName=instrumentName,
2487 refObjectLoader=refObjectLoader,
2488 colorCatalog=colorCatalog,
2489 )
2490 except pipeBase.AlgorithmError as e:
2491 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2492 # No partial outputs for butler to put
2493 raise error from e
2495 for outputRef in outputRefs.outputWcs:
2496 visit = outputRef.dataId["visit"]
2497 butlerQC.put(output.outputWcss[visit], outputRef)
2498 butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
2499 butlerQC.put(output.starCatalog, outputRefs.starCatalog)
2500 if self.config.saveModelParams:
2501 butlerQC.put(output.modelParams, outputRefs.modelParams)
2502 if self.config.saveCameraModel:
2503 butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel)
2504 if self.config.saveCameraObject:
2505 butlerQC.put(output.camera, outputRefs.camera)
2506 if self.config.useColor:
2507 butlerQC.put(output.colorParams, outputRefs.dcrCoefficients)
2508 if output.partialOutputs:
2509 e = RuntimeError("Some visits were dropped because data was insufficient to fit model.")
2510 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2511 raise error from e
2513 def run(
2514 self,
2515 inputVisitSummaries,
2516 isolatedStarSources,
2517 isolatedStarCatalogs,
2518 instrumentName="",
2519 refEpoch=None,
2520 refObjectLoader=None,
2521 inputCameraModel=None,
2522 colorCatalog=None,
2523 inputCamera=None,
2524 ):
2525 """Run the WCS fit for a given set of visits
2527 Parameters
2528 ----------
2529 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
2530 List of catalogs with per-detector summary information.
2531 isolatedStarSources : `list` [`DeferredDatasetHandle`]
2532 List of handles pointing to isolated star sources.
2533 isolatedStarCatalog: `list` [`DeferredDatasetHandle`]
2534 List of handles pointing to isolated star catalogs.
2535 instrumentName : `str`, optional
2536 Name of the instrument used. This is only used for labelling.
2537 refEpoch : `float`, optional
2538 Epoch of the reference objects in MJD.
2539 refObjectLoader : instance of
2540 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`,
2541 optional
2542 Reference object loader instance.
2543 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
2544 Parameters to use for the device part of the model.
2545 colorCatalog : `lsst.afw.table.SimpleCatalog`
2546 Catalog containing object coordinates and magnitudes.
2547 inputCamera : `lsst.afw.cameraGeom.Camera`, optional
2548 Camera to be used as template when constructing new camera.
2550 Returns
2551 -------
2552 result : `lsst.pipe.base.Struct`
2553 ``outputWcss`` : `list` [`lsst.afw.table.ExposureCatalog`]
2554 List of exposure catalogs (one per visit) with the WCS for each
2555 detector set by the new fitted WCS.
2556 ``fitModel`` : `wcsfit.WCSFit`
2557 Model-fitting object with final model parameters.
2558 ``outputCatalog`` : `pyarrow.Table`
2559 Catalog with fit residuals of all sources used.
2560 ``starCatalog`` : `pyarrow.Table`
2561 Catalog with best-fit positions of the objects fit.
2562 ``modelParams`` : `dict`
2563 Parameters and covariance of the best-fit WCS model.
2564 ``cameraModelParams`` : `dict` [`str`, `np.ndarray`]
2565 Parameters of the device part of the model, in the format
2566 needed as input for future runs.
2567 ``colorParams`` : `dict` [`int`, `np.ndarray`]
2568 DCR parameters fit in RA and Dec directions for each visit.
2569 ``camera`` : `lsst.afw.cameraGeom.Camera`
2570 Camera object constructed from the per-detector model.
2571 """
2572 self.log.info("Gather instrument, exposure, and field info")
2574 # Get information about the extent of the input visits
2575 fields, fieldRegions = self._prep_sky(inputVisitSummaries)
2577 # Get RA, Dec, MJD, etc., for the input visits
2578 exposureInfo, exposuresHelper, extensionInfo, instruments = self._get_exposure_info(
2579 inputVisitSummaries, fieldRegions=fieldRegions
2580 )
2582 self.log.info("Load associated sources")
2583 medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="jyear").mjd
2584 allRefObjects, allRefCovariances = {}, {}
2585 for f, fieldRegion in fieldRegions.items():
2586 refObjects, refCovariance = self._load_refcat(
2587 refObjectLoader, extensionInfo, epoch=medianEpoch, region=fieldRegion
2588 )
2589 allRefObjects[f] = refObjects
2590 allRefCovariances[f] = refCovariance
2592 associations, sourceDict = self._associate_from_isolated_sources(
2593 isolatedStarSources, isolatedStarCatalogs, extensionInfo, allRefObjects
2594 )
2596 self.log.info("Fit the WCSs")
2597 # Set up a YAML-type string using the config variables and a sample
2598 # visit
2599 inputYaml, mapTemplate = self.make_yaml(
2600 inputVisitSummaries[0],
2601 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
2602 )
2604 # Set the verbosity level for WCSFit from the task log level.
2605 # TODO: DM-36850, Add lsst.log to gbdes so that log messages are
2606 # properly propagated.
2607 loglevel = self.log.getEffectiveLevel()
2608 if loglevel >= self.log.WARNING:
2609 verbose = 0
2610 elif loglevel == self.log.INFO:
2611 verbose = 1
2612 else:
2613 verbose = 2
2615 # Set up the WCS-fitting class using the source matches from the
2616 # isolated star sources plus the reference catalog.
2617 wcsf = wcsfit.WCSFit(
2618 fields,
2619 instruments,
2620 exposuresHelper,
2621 extensionInfo.visitIndex,
2622 extensionInfo.detectorIndex,
2623 inputYaml,
2624 extensionInfo.wcs,
2625 associations.sequence,
2626 associations.extn,
2627 associations.obj,
2628 sysErr=self.config.systematicError,
2629 refSysErr=self.config.referenceSystematicError,
2630 usePM=self.config.fitProperMotion,
2631 verbose=verbose,
2632 )
2634 # Add the science and reference sources
2635 self._add_objects(wcsf, sourceDict, extensionInfo)
2636 for f in fieldRegions.keys():
2637 self._add_ref_objects(
2638 wcsf, allRefObjects[f], allRefCovariances[f], extensionInfo, fieldIndex=-1 * f
2639 )
2640 if self.config.useColor:
2641 self._add_color_objects(wcsf, colorCatalog)
2643 # Do the WCS fit
2644 try:
2645 wcsf.fit(
2646 reserveFraction=self.config.fitReserveFraction,
2647 randomNumberSeed=self.config.fitReserveRandomSeed,
2648 clipThresh=self.config.clipThresh,
2649 clipFraction=self.config.clipFraction,
2650 )
2651 except RuntimeError as e:
2652 if "Cholesky decomposition failed" in str(e):
2653 raise CholeskyError() from e
2654 else:
2655 raise
2657 self.log.info("WCS fitting done")
2659 outputWcss, cameraParams, colorParams, camera, partialOutputs = self._make_outputs(
2660 wcsf,
2661 inputVisitSummaries,
2662 exposureInfo,
2663 mapTemplate,
2664 extensionInfo,
2665 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
2666 inputCamera=(inputCamera if self.config.buildCamera else None),
2667 )
2668 outputCatalog = wcsf.getOutputCatalog()
2669 outputCatalog["exposureName"] = np.array(outputCatalog["exposureName"])
2670 outputCatalog["deviceName"] = np.array(outputCatalog["deviceName"])
2672 starCatalog = wcsf.getStarCatalog()
2673 modelParams = self._compute_model_params(wcsf) if self.config.saveModelParams else None
2675 return pipeBase.Struct(
2676 outputWcss=outputWcss,
2677 fitModel=wcsf,
2678 outputCatalog=outputCatalog,
2679 starCatalog=starCatalog,
2680 modelParams=modelParams,
2681 cameraModelParams=cameraParams,
2682 colorParams=colorParams,
2683 camera=camera,
2684 partialOutputs=partialOutputs,
2685 )
2687 def _prep_sky(self, inputVisitSummaries):
2688 """Cluster the input visits into disjoint groups that will define
2689 separate fields in the astrometric fit, and, for each group, get the
2690 convex hull around all of its component visits.
2692 The groups are created such that each visit overlaps with at least one
2693 other visit in the same group by the `visitOverlap` amount, which is
2694 calculated as a fraction of the field-of-view radius, and no visits
2695 from separate groups overlap by more than this amount.
2697 Paramaters
2698 ----------
2699 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
2700 List of catalogs with per-detector summary information.
2702 Returns
2703 -------
2704 fields : `wcsfit.Fields`
2705 Object with field information.
2706 fieldRegions : `dict` [`int`, `lsst.sphgeom.ConvexPolygon`]
2707 Dictionary of regions encompassing each group of input visits,
2708 keyed by an arbitrary index.
2709 """
2710 allDetectorCorners = []
2711 mjds = []
2712 radecs = []
2713 radii = []
2714 for visSum in inputVisitSummaries:
2715 detectorCorners = [
2716 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
2717 for (ra, dec) in zip(visSum["raCorners"].ravel(), visSum["decCorners"].ravel())
2718 if (np.isfinite(ra) and (np.isfinite(dec)))
2719 ]
2720 if len(detectorCorners) == 0:
2721 # Skip this visit if none of the detectors have finite ra/dec
2722 # corners, which happens when the WCSs are missing. The visit
2723 # will get formally dropped elsewhere.
2724 continue
2725 allDetectorCorners.append(detectorCorners)
2727 # Get center and approximate radius of field of view
2728 boundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(detectorCorners).getBoundingCircle()
2729 center = lsst.geom.SpherePoint(boundingCircle.getCenter())
2730 ra = center.getRa().asDegrees()
2731 dec = center.getDec().asDegrees()
2732 radecs.append([ra, dec])
2733 radius = boundingCircle.getOpeningAngle()
2734 radii.append(radius)
2736 obsDate = visSum[0].getVisitInfo().getDate()
2737 obsMJD = obsDate.get(obsDate.MJD)
2738 mjds.append(obsMJD)
2740 # Find groups of visits where any one of the visits overlaps another by
2741 # a given fraction of the field-of-view radius.
2742 distance = self.config.visitOverlap * np.median(radii)
2743 clustering = AgglomerativeClustering(
2744 distance_threshold=distance.asDegrees(), n_clusters=None, linkage="single"
2745 )
2746 clusters = clustering.fit(np.array(radecs))
2748 medianMJD = np.median(mjds)
2749 medianEpoch = astropy.time.Time(medianMJD, format="mjd").jyear
2751 fieldNames = []
2752 fieldRAs = []
2753 fieldDecs = []
2754 epochs = []
2755 fieldRegions = {}
2757 for i in range(clusters.n_clusters_):
2758 fieldInd = clusters.labels_ == i
2759 # Concatenate the lists of all detector corners that are in this
2760 # field
2761 fieldDetectors = sum([allDetectorCorners[f] for f, fInd in enumerate(fieldInd) if fInd], [])
2762 hull = lsst.sphgeom.ConvexPolygon.convexHull(fieldDetectors)
2763 center = lsst.geom.SpherePoint(hull.getCentroid())
2764 ra = center.getRa().asDegrees()
2765 dec = center.getDec().asDegrees()
2767 fieldRegions[i] = hull
2768 fieldNames.append(str(i))
2769 fieldRAs.append(ra)
2770 fieldDecs.append(dec)
2771 # Use the same median epoch for all fields so that the final object
2772 # positions are calculated for the same epoch.
2773 epochs.append(medianEpoch)
2775 fields = wcsfit.Fields(fieldNames, fieldRAs, fieldDecs, epochs)
2777 return fields, fieldRegions
2779 def _associate_from_isolated_sources(
2780 self, isolatedStarSourceRefs, isolatedStarCatalogRefs, extensionInfo, refObjects
2781 ):
2782 """Match the input catalog of isolated stars with the reference catalog
2783 and transform the combined isolated star sources and reference source
2784 into the format needed for gbdes.
2786 Parameters
2787 ----------
2788 isolatedStarSourceRefs : `list` [`DeferredDatasetHandle`]
2789 List of handles pointing to isolated star sources.
2790 isolatedStarCatalogRefs: `list` [`DeferredDatasetHandle`]
2791 List of handles pointing to isolated star catalogs.
2792 extensionInfo : `lsst.pipe.base.Struct`
2793 Struct containing properties for each extension (visit/detector).
2794 ``visit`` : `np.ndarray`
2795 Name of the visit for this extension.
2796 ``detector`` : `np.ndarray`
2797 Name of the detector for this extension.
2798 ``visitIndex` : `np.ndarray` [`int`]
2799 Index of visit for this extension.
2800 ``detectorIndex`` : `np.ndarray` [`int`]
2801 Index of the detector for this extension.
2802 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2803 Initial WCS for this extension.
2804 ``extensionType`` : `np.ndarray` [`str`]
2805 "SCIENCE" or "REFERENCE".
2806 refObjects : `dict`
2807 Dictionary of dictionaries containing the position and error
2808 information of reference objects.
2810 Returns
2811 -------
2812 associations : `lsst.pipe.base.Struct`
2813 Struct containing the associations of sources with objects.
2814 sourceDict : `dict` [`int`, [`int`, [`str`, `list` [`float`]]]]
2815 Dictionary containing the source centroids for each visit.
2816 """
2817 sequences = []
2818 extensions = []
2819 object_indices = []
2821 sourceColumns = ["x", "y", "xErr", "yErr", "ixx", "ixy", "iyy", "obj_index", "visit", "detector"]
2822 catalogColumns = ["ra", "dec"]
2824 sourceDict = dict([(visit, {}) for visit in np.unique(extensionInfo.visit)])
2825 for visit, detector in zip(extensionInfo.visit, extensionInfo.detector):
2826 sourceDict[visit][detector] = {"x": [], "y": [], "xCov": [], "yCov": [], "xyCov": []}
2828 for isolatedStarCatalogRef, isolatedStarSourceRef in zip(
2829 isolatedStarCatalogRefs, isolatedStarSourceRefs
2830 ):
2831 isolatedStarCatalog = isolatedStarCatalogRef.get(parameters={"columns": catalogColumns})
2832 isolatedStarSources = isolatedStarSourceRef.get(parameters={"columns": sourceColumns})
2833 if len(isolatedStarCatalog) == 0:
2834 # This is expected when only one visit overlaps with a given
2835 # tract, meaning that no sources can be associated.
2836 self.log.debug(
2837 "Skipping tract %d, which has no associated isolated stars",
2838 isolatedStarCatalogRef.dataId["tract"],
2839 )
2840 continue
2842 # Match the reference stars to the existing isolated stars, then
2843 # insert the reference stars into the isolated star sources.
2844 allVisits = np.copy(isolatedStarSources["visit"])
2845 allDetectors = np.copy(isolatedStarSources["detector"])
2846 allObjectIndices = np.copy(isolatedStarSources["obj_index"])
2847 issIndices = np.copy(isolatedStarSources.index)
2848 for f, regionRefObjects in refObjects.items():
2849 # Use the same matching technique that is done in
2850 # isolatedStarAssociation and fgcmBuildFromIsolatedStars.
2851 with Matcher(
2852 isolatedStarCatalog["ra"].to_numpy(), isolatedStarCatalog["dec"].to_numpy()
2853 ) as matcher:
2854 idx, idx_isoStarCat, idx_refObjects, d = matcher.query_radius(
2855 np.array(regionRefObjects["ra"]),
2856 np.array(regionRefObjects["dec"]),
2857 self.config.matchRadius / 3600.0,
2858 return_indices=True,
2859 )
2861 refSort = np.searchsorted(isolatedStarSources["obj_index"], idx_isoStarCat)
2862 refDetector = np.ones(len(idx_isoStarCat)) * -1
2863 # The "visit" for the reference catalogs is the field times -1.
2864 refVisit = np.ones(len(idx_isoStarCat)) * f * -1
2866 allVisits = np.insert(allVisits, refSort, refVisit)
2867 allDetectors = np.insert(allDetectors, refSort, refDetector)
2868 allObjectIndices = np.insert(allObjectIndices, refSort, idx_isoStarCat)
2869 issIndices = np.insert(issIndices, refSort, idx_refObjects)
2871 # Loop through the associated sources to convert them to the gbdes
2872 # format, which requires the extension index, the source's index in
2873 # the input table, and a sequence number corresponding to the
2874 # object with which it is associated.
2875 sequence = 0
2876 obj_index = allObjectIndices[0]
2877 for visit, detector, row, obj_ind in zip(allVisits, allDetectors, issIndices, allObjectIndices):
2878 extensionIndex = np.flatnonzero(
2879 (extensionInfo.visit == visit) & (extensionInfo.detector == detector)
2880 )
2881 if len(extensionIndex) == 0:
2882 # This happens for runs where you are not using all the
2883 # visits overlapping a given tract that were included in
2884 # the isolated star association task."
2885 continue
2886 else:
2887 extensionIndex = extensionIndex[0]
2889 extensions.append(extensionIndex)
2890 if visit <= 0:
2891 object_indices.append(row)
2892 else:
2893 object_indices.append(len(sourceDict[visit][detector]["x"]))
2894 source = isolatedStarSources.loc[row]
2895 sourceDict[visit][detector]["x"].append(source["x"])
2896 sourceDict[visit][detector]["y"].append(source["y"])
2897 xCov = source["xErr"] ** 2
2898 yCov = source["yErr"] ** 2
2899 xyCov = source["ixy"] * (xCov + yCov) / (source["ixx"] + source["iyy"])
2900 # TODO: add correct xyErr if DM-7101 is ever done.
2901 sourceDict[visit][detector]["xCov"].append(xCov)
2902 sourceDict[visit][detector]["yCov"].append(yCov)
2903 sourceDict[visit][detector]["xyCov"].append(xyCov)
2904 if obj_ind != obj_index:
2905 sequence = 0
2906 sequences.append(sequence)
2907 obj_index = obj_ind
2908 sequence += 1
2909 else:
2910 sequences.append(sequence)
2911 sequence += 1
2913 associations = pipeBase.Struct(extn=extensions, obj=object_indices, sequence=sequences)
2914 return associations, sourceDict
2916 def _add_objects(self, wcsf, sourceDict, extensionInfo):
2917 """Add science sources to the wcsfit.WCSFit object.
2919 Parameters
2920 ----------
2921 wcsf : `wcsfit.WCSFit`
2922 WCS-fitting object.
2923 sourceDict : `dict`
2924 Dictionary containing the source centroids for each visit.
2925 extensionInfo : `lsst.pipe.base.Struct`
2926 Struct containing properties for each extension (visit/detector).
2927 ``visit`` : `np.ndarray`
2928 Name of the visit for this extension.
2929 ``detector`` : `np.ndarray`
2930 Name of the detector for this extension.
2931 ``visitIndex` : `np.ndarray` [`int`]
2932 Index of visit for this extension.
2933 ``detectorIndex`` : `np.ndarray` [`int`]
2934 Index of the detector for this extension.
2935 ``wcss`` : `np.ndarray` [`lsst.afw.geom.SkyWcs`]
2936 Initial WCS for this extension.
2937 ``extensionType`` : `np.ndarray` [`str`]
2938 "SCIENCE" or "REFERENCE".
2939 """
2940 for visit, visitSources in sourceDict.items():
2941 # Visit numbers equal or below zero connote the reference catalog.
2942 if visit <= 0:
2943 # This "visit" number corresponds to a reference catalog.
2944 continue
2946 for detector, sourceCat in visitSources.items():
2947 extensionIndex = np.flatnonzero(
2948 (extensionInfo.visit == visit) & (extensionInfo.detector == detector)
2949 )[0]
2951 d = {
2952 "x": np.array(sourceCat["x"]),
2953 "y": np.array(sourceCat["y"]),
2954 "xCov": np.array(sourceCat["xCov"]),
2955 "yCov": np.array(sourceCat["yCov"]),
2956 "xyCov": np.array(sourceCat["xyCov"]),
2957 }
2958 wcsf.setObjects(
2959 extensionIndex,
2960 d,
2961 "x",
2962 "y",
2963 ["xCov", "yCov", "xyCov"],
2964 defaultColor=self.config.referenceColor,
2965 )
2968class GbdesGlobalAstrometricMultibandFitConnections(
2969 GbdesGlobalAstrometricFitConnections,
2970 dimensions=("instrument",),
2971):
2972 outputCatalog = pipeBase.connectionTypes.Output(
2973 doc=(
2974 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2975 "plane coordinates and chisq values."
2976 ),
2977 name="gbdesGlobalAstrometricMultibandFit_fitStars",
2978 storageClass="ArrowNumpyDict",
2979 dimensions=("instrument",),
2980 )
2981 starCatalog = pipeBase.connectionTypes.Output(
2982 doc=(
2983 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2984 "fitProperMotion is True."
2985 ),
2986 name="gbdesGlobalAstrometricMultibandFit_starCatalog",
2987 storageClass="ArrowNumpyDict",
2988 dimensions=("instrument",),
2989 )
2990 modelParams = pipeBase.connectionTypes.Output(
2991 doc="WCS parameters and covariance.",
2992 name="gbdesGlobalAstrometricMultibandFit_modelParams",
2993 storageClass="ArrowNumpyDict",
2994 dimensions=("instrument",),
2995 )
2998class GbdesGlobalAstrometricMultibandFitConfig(
2999 GbdesAstrometricFitConfig,
3000 pipelineConnections=GbdesGlobalAstrometricMultibandFitConnections,
3001):
3002 """Configuration for the GbdesGlobalAstrometricMultibandFitTask"""
3004 pass
3007class GbdesGlobalAstrometricMultibandFitTask(GbdesGlobalAstrometricFitTask):
3008 """Calibrate the WCS across multiple visits in multiple filters and
3009 multiple fields using the GBDES package.
3010 """
3012 ConfigClass = GbdesGlobalAstrometricMultibandFitConfig
3013 _DefaultName = "gbdesAstrometricMultibandFit"