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