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