Coverage for python/lsst/drp/tasks/gbdesAstrometricFit.py: 71%
1111 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-15 09:58 +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"]: 790 ↛ 791line 790 didn't jump to line 791 because the condition on line 790 was never true
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"]: 799 ↛ 800line 799 didn't jump to line 800 because the condition on line 799 was never true
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: 806 ↛ 807line 806 didn't jump to line 807 because the condition on line 806 was never true
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): 813 ↛ 814line 813 didn't jump to line 814 because the condition on line 813 was never true
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: 833 ↛ 834line 833 didn't jump to line 834 because the condition on line 833 was never true
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: 1036 ↛ 1037line 1036 didn't jump to line 1037 because the condition on line 1036 was never true
1037 verbose = 0
1038 elif loglevel == self.log.INFO: 1038 ↛ 1039line 1038 didn't jump to line 1039 because the condition on line 1038 was never true
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: 1255 ↛ 1256line 1255 didn't jump to line 1256 because the condition on line 1255 was never true
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: 1280 ↛ 1281line 1280 didn't jump to line 1281 because the condition on line 1280 was never true
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: 1324 ↛ 1325line 1324 didn't jump to line 1325 because the condition on line 1324 was never true
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: 1329 ↛ 1333line 1329 didn't jump to line 1333 because the condition on line 1329 was always true
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: 1347 ↛ 1348line 1347 didn't jump to line 1348 because the condition on line 1347 was never true
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: 1452 ↛ 1453line 1452 didn't jump to line 1453 because the condition on line 1452 was never true
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: 1458 ↛ 1459line 1458 didn't jump to line 1459 because the condition on line 1458 was never true
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: 1471 ↛ 1472line 1471 didn't jump to line 1472 because the condition on line 1471 was never true
1472 refCatVersion = 2
1473 else:
1474 refCatVersion = refCatMetadata["REFCAT_FORMAT_VERSION"]
1475 if refCatVersion == 2: 1475 ↛ 1476line 1475 didn't jump to line 1476 because the condition on line 1475 was never true
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): 1504 ↛ 1509line 1504 didn't jump to line 1509 because the condition on line 1504 was always true
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(): 1515 ↛ 1516line 1515 didn't jump to line 1516 because the condition on line 1515 was never true
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: 1527 ↛ 1529line 1527 didn't jump to line 1529 because the condition on line 1527 was never true
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: 1539 ↛ 1540line 1539 didn't jump to line 1540 because the condition on line 1539 was never true
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: 1547 ↛ 1548line 1547 didn't jump to line 1548 because the condition on line 1547 was never true
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: 1662 ↛ 1664line 1662 didn't jump to line 1664 because the condition on line 1662 was always true
1663 columns.extend(self.sourceSelector.config.flags.bad)
1664 if self.sourceSelector.config.doUnresolved: 1664 ↛ 1666line 1664 didn't jump to line 1666 because the condition on line 1664 was always true
1665 columns.append(self.sourceSelector.config.unresolved.name)
1666 if self.sourceSelector.config.doIsolated: 1666 ↛ 1669line 1666 didn't jump to line 1669 because the condition on line 1666 was always true
1667 columns.append(self.sourceSelector.config.isolated.parentName)
1668 columns.append(self.sourceSelector.config.isolated.nChildName)
1669 if self.sourceSelector.config.doRequirePrimary: 1669 ↛ 1672line 1669 didn't jump to line 1672 because the condition on line 1669 was always true
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: 1781 ↛ 1784line 1781 didn't jump to line 1784 because the condition on line 1781 was never true
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(): 1797 ↛ 1798line 1797 didn't jump to line 1798 because the condition on line 1797 was never true
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)): 1822 ↛ 1823line 1822 didn't jump to line 1823 because the condition on line 1822 was never true
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: 1953 ↛ 1954line 1953 didn't jump to line 1954 because the condition on line 1953 was never true
1954 unconstrainedDetectors.append(str(detector))
1956 if unconstrainedDetectors: 1956 ↛ 1957line 1956 didn't jump to line 1957 because the condition on line 1956 was never true
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: 1983 ↛ 1984line 1983 didn't jump to line 1984 because the condition on line 1983 was never true
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": 2001 ↛ 2011line 2001 didn't jump to line 2011 because the condition on line 2001 was always true
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: 2020 ↛ 2021line 2020 didn't jump to line 2021 because the condition on line 2020 was never true
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": 2056 ↛ 2057line 2056 didn't jump to line 2057 because the condition on line 2056 was never true
2057 componentDict = {"Type": "Identity"}
2058 elif mapping == "dcr": 2058 ↛ 2065line 2058 didn't jump to line 2065 because the condition on line 2058 was always true
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: 2222 ↛ 2223line 2222 didn't jump to line 2223 because the condition on line 2222 was never true
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: 2304 ↛ 2309line 2304 didn't jump to line 2309 because the condition on line 2304 was always true
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": 2330 ↛ 2333line 2330 didn't jump to line 2333 because the condition on line 2330 was always true
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: 2412 ↛ 2414line 2412 didn't jump to line 2414 because the condition on line 2412 was never true
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: 2433 ↛ 2434line 2433 didn't jump to line 2434 because the condition on line 2433 was never true
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: 2455 ↛ 2454line 2455 didn't jump to line 2454 because the loop on line 2455 didn't complete
2456 if detInfo := visitSummary.find(detector): 2456 ↛ 2455line 2456 didn't jump to line 2455 because the condition on line 2456 was always true
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: 2465 ↛ 2466line 2465 didn't jump to line 2466 because the condition on line 2465 was never true
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 visitMap = visitMaps[0]
2472 visitMapType = wcsf.mapCollection.getMapType(visitMap)
2473 if (visitMap not in mapParams) and (visitMapType != "Identity"): 2473 ↛ 2474line 2473 didn't jump to line 2474 because the condition on line 2473 was never true
2474 self.log.warning("Visit %d was dropped because of an insufficient amount of data.", visit)
2475 partialOutputs = True
2476 continue
2478 if self.config.useColor:
2479 if f"{visit}/dcr" not in visitMaps: 2479 ↛ 2480line 2479 didn't jump to line 2480 because the condition on line 2479 was never true
2480 raise RuntimeError("DCR map missing for visit %s", visit)
2481 colorMap = visitMaps.pop(visitMaps.index(f"{visit}/dcr"))
2482 colorFits[visit] = mapParams[colorMap]
2484 catalog = lsst.afw.table.ExposureCatalog(schema)
2485 catalog.resize(len(exposureInfo.detectors))
2486 catalog["visit"] = visit
2488 for d, detector in enumerate(exposureInfo.detectors):
2489 mapName = f"{visit}/{detector}"
2490 if mapName in wcsf.mapCollection.allMapNames():
2491 mapElements = wcsf.mapCollection.orderAtoms(f"{mapName}/base")
2492 catalog[d]["recoveredWcs"] = False
2493 else:
2494 # This extension was not fit, but the WCS can be recovered
2495 # using the maps fit from sources on other visits but the
2496 # same detector and from sources on other detectors from
2497 # this visit.
2498 genericElements = mapTemplate["EXPOSURE/DEVICE/base"]["Elements"]
2499 mapElements = []
2500 band = visitSummary[0]["physical_filter"]
2502 # Go through the generic map components to build the names
2503 # of the specific maps for this extension.
2504 for component in genericElements:
2505 elements = mapTemplate[component]["Elements"]
2506 for element in elements:
2507 # TODO: DM-42519, gbdes sets the "BAND" to the
2508 # instrument name currently. This will need to be
2509 # disambiguated if we run on multiple bands at
2510 # once.
2511 element = element.replace("BAND", str(band))
2512 element = element.replace("EXPOSURE", str(visit))
2513 element = element.replace("DEVICE", str(detector))
2514 mapElements.append(element)
2515 catalog[d]["recoveredWcs"] = True
2516 mapDict = {}
2517 for m, mapElement in enumerate(mapElements):
2518 mapType = wcsf.mapCollection.getMapType(mapElement)
2519 if mapType == "Color":
2520 # DCR fit should not go into the generic WCS.
2521 continue
2522 mapDict[mapElement] = {"Type": mapType}
2524 if mapType == "Poly":
2525 mapCoefficients = mapParams[mapElement]
2526 mapDict[mapElement]["Coefficients"] = mapCoefficients
2528 # The RA and Dec of the visit are needed for the last step of
2529 # the mapping from the visit tangent plane to RA and Dec
2530 outWCS = self._make_afw_wcs(
2531 mapDict,
2532 exposureInfo.ras[v] * lsst.geom.radians,
2533 exposureInfo.decs[v] * lsst.geom.radians,
2534 doNormalizePixels=True,
2535 xScale=xscale,
2536 yScale=yscale,
2537 )
2539 catalog[d].setId(detector)
2540 catalog[d].setBBox(detectorBBoxes[detector])
2541 catalog[d].setWcs(outWCS)
2542 catalog.sort()
2543 catalogs[visit] = catalog
2544 if self.config.useColor:
2545 colorVisits = np.array(list(colorFits.keys()))
2546 colorRA = np.array([colorFits[vis][0] for vis in colorVisits])
2547 colorDec = np.array([colorFits[vis][1] for vis in colorVisits])
2548 colorFits = {"visit": colorVisits, "raCoefficient": colorRA, "decCoefficient": colorDec}
2550 return catalogs, cameraParams, colorFits, camera, partialOutputs
2552 def _compute_model_params(self, wcsf):
2553 """Get the WCS model parameters and covariance and convert to a
2554 dictionary that will be readable as a pandas dataframe or other table.
2556 Parameters
2557 ----------
2558 wcsf : `wcsfit.WCSFit`
2559 WCSFit object, assumed to have fit model.
2561 Returns
2562 -------
2563 modelParams : `dict`
2564 Parameters and covariance of the best-fit WCS model.
2565 """
2566 modelParamDict = wcsf.mapCollection.getParamDict()
2567 modelCovariance = wcsf.getModelCovariance()
2569 modelParams = {k: [] for k in ["mapName", "coordinate", "parameter", "coefficientNumber"]}
2570 i = 0
2571 for mapName, params in modelParamDict.items():
2572 nCoeffs = len(params)
2573 # There are an equal number of x and y coordinate parameters
2574 nCoordCoeffs = nCoeffs // 2
2575 modelParams["mapName"].extend([mapName] * nCoeffs)
2576 modelParams["coordinate"].extend(["x"] * nCoordCoeffs)
2577 modelParams["coordinate"].extend(["y"] * nCoordCoeffs)
2578 modelParams["parameter"].extend(params)
2579 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2580 modelParams["coefficientNumber"].extend(np.arange(nCoordCoeffs))
2582 for p in range(nCoeffs):
2583 if p < nCoordCoeffs:
2584 coord = "x"
2585 else:
2586 coord = "y"
2587 modelParams[f"{mapName}_{coord}_{p}_cov"] = modelCovariance[i]
2588 i += 1
2590 # Convert the dictionary values from lists to numpy arrays.
2591 for key, value in modelParams.items():
2592 modelParams[key] = np.array(value)
2594 return modelParams
2597class GbdesAstrometricMultibandFitConnections(
2598 GbdesAstrometricFitConnections, dimensions=("skymap", "tract", "instrument")
2599):
2600 outputCatalog = pipeBase.connectionTypes.Output(
2601 doc=(
2602 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2603 "plane coordinates and chisq values."
2604 ),
2605 name="gbdesAstrometricMultibandFit_fitStars",
2606 storageClass="ArrowNumpyDict",
2607 dimensions=("instrument", "skymap", "tract"),
2608 )
2609 starCatalog = pipeBase.connectionTypes.Output(
2610 doc=(
2611 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2612 "fitProperMotion is True."
2613 ),
2614 name="gbdesAstrometricMultibandFit_starCatalog",
2615 storageClass="ArrowNumpyDict",
2616 dimensions=("instrument", "skymap", "tract"),
2617 )
2618 modelParams = pipeBase.connectionTypes.Output(
2619 doc="WCS parameters and covariance.",
2620 name="gbdesAstrometricMultibandFit_modelParams",
2621 storageClass="ArrowNumpyDict",
2622 dimensions=("instrument", "skymap", "tract"),
2623 )
2626class GbdesAstrometricMultibandFitConfig(
2627 GbdesAstrometricFitConfig, pipelineConnections=GbdesAstrometricMultibandFitConnections
2628):
2629 pass
2632class GbdesAstrometricMultibandFitTask(GbdesAstrometricFitTask):
2633 """Calibrate the WCS across multiple visits in multiple filters of the same
2634 field using the GBDES package.
2635 """
2637 ConfigClass = GbdesAstrometricMultibandFitConfig
2638 _DefaultName = "gbdesAstrometricMultibandFit"
2641class GbdesGlobalAstrometricFitConnections(
2642 pipeBase.PipelineTaskConnections, dimensions=("instrument", "physical_filter")
2643):
2644 inputVisitSummaries = pipeBase.connectionTypes.Input(
2645 doc=(
2646 "Per-visit consolidated exposure metadata built from calexps. "
2647 "These catalogs use detector id for the id and must be sorted for "
2648 "fast lookups of a detector."
2649 ),
2650 name="visitSummary",
2651 storageClass="ExposureCatalog",
2652 dimensions=("instrument", "visit"),
2653 multiple=True,
2654 )
2655 referenceCatalog = pipeBase.connectionTypes.PrerequisiteInput(
2656 doc="The astrometry reference catalog to match to loaded input catalog sources.",
2657 name="gaia_dr3_20230707",
2658 storageClass="SimpleCatalog",
2659 dimensions=("skypix",),
2660 deferLoad=True,
2661 multiple=True,
2662 )
2663 colorCatalog = pipeBase.connectionTypes.Input(
2664 doc="The catalog of magnitudes to match to input sources.",
2665 name="fgcm_Cycle4_StandardStars",
2666 storageClass="SimpleCatalog",
2667 dimensions=("instrument",),
2668 )
2669 isolatedStarSources = pipeBase.connectionTypes.Input(
2670 doc="Catalog of matched sources.",
2671 name="isolated_star_presources",
2672 storageClass="DataFrame",
2673 dimensions=(
2674 "instrument",
2675 "skymap",
2676 "tract",
2677 ),
2678 multiple=True,
2679 deferLoad=True,
2680 )
2681 isolatedStarCatalogs = pipeBase.connectionTypes.Input(
2682 doc="Catalog of objects corresponding to the isolatedStarSources.",
2683 name="isolated_star_presource_associations",
2684 storageClass="DataFrame",
2685 dimensions=(
2686 "instrument",
2687 "skymap",
2688 "tract",
2689 ),
2690 multiple=True,
2691 deferLoad=True,
2692 )
2693 inputCameraModel = pipeBase.connectionTypes.PrerequisiteInput(
2694 doc="Camera parameters to use for 'device' part of model",
2695 name="gbdesAstrometricFit_cameraModel",
2696 storageClass="ArrowNumpyDict",
2697 dimensions=("instrument", "physical_filter"),
2698 )
2699 inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
2700 doc="Input camera to use when constructing camera from astrometric model.",
2701 name="camera",
2702 storageClass="Camera",
2703 dimensions=("instrument",),
2704 isCalibration=True,
2705 )
2706 outputWcs = pipeBase.connectionTypes.Output(
2707 doc=(
2708 "Per-visit world coordinate systems derived from the fitted model. These catalogs only contain "
2709 "entries for detectors with an output, and use the detector id for the catalog id, sorted on id "
2710 "for fast lookups of a detector."
2711 ),
2712 name="gbdesGlobalAstrometricFitSkyWcsCatalog",
2713 storageClass="ExposureCatalog",
2714 dimensions=("instrument", "visit"),
2715 multiple=True,
2716 )
2717 outputCatalog = pipeBase.connectionTypes.Output(
2718 doc=(
2719 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
2720 "plane coordinates and chisq values."
2721 ),
2722 name="gbdesGlobalAstrometricFit_fitStars",
2723 storageClass="ArrowNumpyDict",
2724 dimensions=("instrument", "physical_filter"),
2725 )
2726 starCatalog = pipeBase.connectionTypes.Output(
2727 doc=(
2728 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
2729 "fitProperMotion is True."
2730 ),
2731 name="gbdesGlobalAstrometricFit_starCatalog",
2732 storageClass="ArrowNumpyDict",
2733 dimensions=("instrument", "physical_filter"),
2734 )
2735 modelParams = pipeBase.connectionTypes.Output(
2736 doc="WCS parameters and covariance.",
2737 name="gbdesGlobalAstrometricFit_modelParams",
2738 storageClass="ArrowNumpyDict",
2739 dimensions=("instrument", "physical_filter"),
2740 )
2741 outputCameraModel = pipeBase.connectionTypes.Output(
2742 doc="Camera parameters to use for 'device' part of model",
2743 name="gbdesAstrometricFit_cameraModel",
2744 storageClass="ArrowNumpyDict",
2745 dimensions=("instrument", "physical_filter"),
2746 )
2747 dcrCoefficients = pipeBase.connectionTypes.Output(
2748 doc="Per-visit coefficients for DCR correction.",
2749 name="gbdesGlobalAstrometricFit_dcrCoefficients",
2750 storageClass="ArrowNumpyDict",
2751 )
2752 camera = pipeBase.connectionTypes.Output(
2753 doc="Camera object constructed using the per-detector part of the astrometric model",
2754 name="gbdesGlobalAstrometricFitCamera",
2755 storageClass="Camera",
2756 dimensions=("instrument", "physical_filter"),
2757 )
2759 def getSpatialBoundsConnections(self):
2760 return ("inputVisitSummaries",)
2762 def __init__(self, *, config=None):
2763 super().__init__(config=config)
2765 if not self.config.useColor:
2766 self.inputs.remove("colorCatalog")
2767 self.outputs.remove("dcrCoefficients")
2768 if not self.config.saveModelParams:
2769 self.outputs.remove("modelParams")
2770 if not self.config.useInputCameraModel:
2771 self.prerequisiteInputs.remove("inputCameraModel")
2772 if not self.config.saveCameraModel:
2773 self.outputs.remove("outputCameraModel")
2774 if not self.config.saveCameraObject:
2775 self.prerequisiteInputs.remove("inputCamera")
2776 self.outputs.remove("camera")
2779class GbdesGlobalAstrometricFitConfig(
2780 GbdesAstrometricFitConfig, pipelineConnections=GbdesGlobalAstrometricFitConnections
2781):
2782 visitOverlap = pexConfig.Field(
2783 dtype=float,
2784 default=1.0,
2785 doc=(
2786 "The linkage distance threshold above which clustered groups of visits will not be merged "
2787 "together in an agglomerative clustering algorithm. The linkage distance is calculated using the "
2788 "minimum distance between the field-of-view centers of a given visit and all other visits in a "
2789 "group, and is in units of the field-of-view radius. The resulting groups of visits define the "
2790 "fields for the astrometric fit."
2791 ),
2792 )
2795class GbdesGlobalAstrometricFitTask(GbdesAstrometricFitTask):
2796 """Calibrate the WCS across multiple visits and multiple fields using the
2797 GBDES package.
2799 This class assumes that the input visits can be separated into contiguous
2800 groups, for which an individual group covers an area of less than a
2801 hemisphere.
2802 """
2804 ConfigClass = GbdesGlobalAstrometricFitConfig
2805 _DefaultName = "gbdesAstrometricFit"
2807 def runQuantum(self, butlerQC, inputRefs, outputRefs):
2808 # We override runQuantum to set up the refObjLoaders
2809 inputs = butlerQC.get(inputRefs)
2811 instrumentName = butlerQC.quantum.dataId["instrument"]
2813 # Ensure the inputs are in a consistent and deterministic order
2814 inputSumVisits = np.array([inputSum[0]["visit"] for inputSum in inputs["inputVisitSummaries"]])
2815 inputs["inputVisitSummaries"] = [inputs["inputVisitSummaries"][v] for v in inputSumVisits.argsort()]
2816 inputRefHtm7s = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputRefs.referenceCatalog])
2817 inputRefCatRefs = [inputRefs.referenceCatalog[htm7] for htm7 in inputRefHtm7s.argsort()]
2818 inputRefCats = np.array([inputRefCat.dataId["htm7"] for inputRefCat in inputs["referenceCatalog"]])
2819 inputs["referenceCatalog"] = [inputs["referenceCatalog"][v] for v in inputRefCats.argsort()]
2820 inputIsolatedStarSourceTracts = np.array(
2821 [isolatedStarSource.dataId["tract"] for isolatedStarSource in inputs["isolatedStarSources"]]
2822 )
2823 inputIsolatedStarCatalogTracts = np.array(
2824 [isolatedStarCatalog.dataId["tract"] for isolatedStarCatalog in inputs["isolatedStarCatalogs"]]
2825 )
2826 for tract in inputIsolatedStarCatalogTracts:
2827 if tract not in inputIsolatedStarSourceTracts:
2828 raise RuntimeError(f"tract {tract} in isolated_star_cats but not isolated_star_sources")
2829 inputs["isolatedStarSources"] = np.array(
2830 [inputs["isolatedStarSources"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2831 )
2832 inputs["isolatedStarCatalogs"] = np.array(
2833 [inputs["isolatedStarCatalogs"][t] for t in inputIsolatedStarSourceTracts.argsort()]
2834 )
2836 refConfig = LoadReferenceObjectsConfig()
2837 if self.config.applyRefCatProperMotion:
2838 refConfig.requireProperMotion = True
2839 refObjectLoader = ReferenceObjectLoader(
2840 dataIds=[ref.datasetRef.dataId for ref in inputRefCatRefs],
2841 refCats=inputs.pop("referenceCatalog"),
2842 config=refConfig,
2843 log=self.log,
2844 )
2845 if self.config.useColor:
2846 colorCatalog = inputs.pop("colorCatalog")
2847 else:
2848 colorCatalog = None
2850 try:
2851 output = self.run(
2852 **inputs,
2853 instrumentName=instrumentName,
2854 refObjectLoader=refObjectLoader,
2855 colorCatalog=colorCatalog,
2856 )
2857 except pipeBase.AlgorithmError as e:
2858 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2859 # No partial outputs for butler to put
2860 raise error from e
2862 for outputRef in outputRefs.outputWcs:
2863 visit = outputRef.dataId["visit"]
2864 butlerQC.put(output.outputWcss[visit], outputRef)
2865 butlerQC.put(output.outputCatalog, outputRefs.outputCatalog)
2866 butlerQC.put(output.starCatalog, outputRefs.starCatalog)
2867 if self.config.saveModelParams:
2868 butlerQC.put(output.modelParams, outputRefs.modelParams)
2869 if self.config.saveCameraModel:
2870 butlerQC.put(output.cameraModelParams, outputRefs.outputCameraModel)
2871 if self.config.saveCameraObject:
2872 butlerQC.put(output.camera, outputRefs.camera)
2873 if self.config.useColor:
2874 butlerQC.put(output.colorParams, outputRefs.dcrCoefficients)
2875 if output.partialOutputs:
2876 e = RuntimeError("Some visits were dropped because data was insufficient to fit model.")
2877 error = pipeBase.AnnotatedPartialOutputsError.annotate(e, self, log=self.log)
2878 raise error from e
2880 def run(
2881 self,
2882 inputVisitSummaries,
2883 isolatedStarSources,
2884 isolatedStarCatalogs,
2885 instrumentName="",
2886 refEpoch=None,
2887 refObjectLoader=None,
2888 inputCameraModel=None,
2889 colorCatalog=None,
2890 inputCamera=None,
2891 ):
2892 """Run the WCS fit for a given set of visits
2894 Parameters
2895 ----------
2896 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
2897 List of catalogs with per-detector summary information.
2898 isolatedStarSources : `list` [`DeferredDatasetHandle`]
2899 List of handles pointing to isolated star sources.
2900 isolatedStarCatalog: `list` [`DeferredDatasetHandle`]
2901 List of handles pointing to isolated star catalogs.
2902 instrumentName : `str`, optional
2903 Name of the instrument used. This is only used for labelling.
2904 refEpoch : `float`, optional
2905 Epoch of the reference objects in MJD.
2906 refObjectLoader : instance of
2907 `lsst.meas.algorithms.loadReferenceObjects.ReferenceObjectLoader`,
2908 optional
2909 Reference object loader instance.
2910 inputCameraModel : `dict` [`str`, `np.ndarray`], optional
2911 Parameters to use for the device part of the model.
2912 colorCatalog : `lsst.afw.table.SimpleCatalog`
2913 Catalog containing object coordinates and magnitudes.
2914 inputCamera : `lsst.afw.cameraGeom.Camera`, optional
2915 Camera to be used as template when constructing new camera.
2917 Returns
2918 -------
2919 result : `lsst.pipe.base.Struct`
2920 ``outputWcss`` : `list` [`lsst.afw.table.ExposureCatalog`]
2921 List of exposure catalogs (one per visit) with the WCS for each
2922 detector set by the new fitted WCS.
2923 ``fitModel`` : `wcsfit.WCSFit`
2924 Model-fitting object with final model parameters.
2925 ``outputCatalog`` : `pyarrow.Table`
2926 Catalog with fit residuals of all sources used.
2927 ``starCatalog`` : `pyarrow.Table`
2928 Catalog with best-fit positions of the objects fit.
2929 ``modelParams`` : `dict`
2930 Parameters and covariance of the best-fit WCS model.
2931 ``cameraModelParams`` : `dict` [`str`, `np.ndarray`]
2932 Parameters of the device part of the model, in the format
2933 needed as input for future runs.
2934 ``colorParams`` : `dict` [`int`, `np.ndarray`]
2935 DCR parameters fit in RA and Dec directions for each visit.
2936 ``camera`` : `lsst.afw.cameraGeom.Camera`
2937 Camera object constructed from the per-detector model.
2938 """
2939 self.log.info("Gather instrument, exposure, and field info")
2941 # Get information about the extent of the input visits
2942 fields, fieldRegions = self._prep_sky(inputVisitSummaries)
2944 # Get RA, Dec, MJD, etc., for the input visits
2945 exposureInfo, exposuresHelper, extensionInfo, instruments = self._get_exposure_info(
2946 inputVisitSummaries, fieldRegions=fieldRegions
2947 )
2949 self.log.info("Load associated sources")
2950 medianEpoch = astropy.time.Time(exposureInfo.medianEpoch, format="jyear").mjd
2951 allRefObjects, allRefCovariances = {}, {}
2952 for f, fieldRegion in fieldRegions.items():
2953 refObjects, refCovariance = self._load_refcat(
2954 refObjectLoader, extensionInfo, epoch=medianEpoch, region=fieldRegion
2955 )
2956 allRefObjects[f] = refObjects
2957 allRefCovariances[f] = refCovariance
2959 associations, sourceDict = self._associate_from_isolated_sources(
2960 isolatedStarSources, isolatedStarCatalogs, extensionInfo, allRefObjects
2961 )
2963 self.log.info("Fit the WCSs")
2964 # Set up a YAML-type string using the config variables and a sample
2965 # visit
2966 inputYaml, mapTemplate = self.make_yaml(
2967 inputVisitSummaries[0],
2968 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
2969 )
2971 # Set the verbosity level for WCSFit from the task log level.
2972 # TODO: DM-36850, Add lsst.log to gbdes so that log messages are
2973 # properly propagated.
2974 loglevel = self.log.getEffectiveLevel()
2975 if loglevel >= self.log.WARNING: 2975 ↛ 2976line 2975 didn't jump to line 2976 because the condition on line 2975 was never true
2976 verbose = 0
2977 elif loglevel == self.log.INFO: 2977 ↛ 2978line 2977 didn't jump to line 2978 because the condition on line 2977 was never true
2978 verbose = 1
2979 else:
2980 verbose = 2
2982 # Set up the WCS-fitting class using the source matches from the
2983 # isolated star sources plus the reference catalog.
2984 wcsf = wcsfit.WCSFit(
2985 fields,
2986 instruments,
2987 exposuresHelper,
2988 extensionInfo.visitIndex,
2989 extensionInfo.detectorIndex,
2990 inputYaml,
2991 extensionInfo.wcs,
2992 associations.sequence,
2993 associations.extn,
2994 associations.obj,
2995 sysErr=self.config.systematicError,
2996 refSysErr=self.config.referenceSystematicError,
2997 usePM=self.config.fitProperMotion,
2998 verbose=verbose,
2999 )
3001 # Add the science and reference sources
3002 self._add_objects_from_dict(wcsf, sourceDict, extensionInfo)
3003 for f in fieldRegions.keys():
3004 self._add_ref_objects(
3005 wcsf, allRefObjects[f], allRefCovariances[f], extensionInfo, fieldIndex=-1 * f
3006 )
3007 if self.config.useColor:
3008 self._add_color_objects(wcsf, colorCatalog)
3010 # Do the WCS fit
3011 try:
3012 wcsf.fit(
3013 reserveFraction=self.config.fitReserveFraction,
3014 randomNumberSeed=self.config.fitReserveRandomSeed,
3015 clipThresh=self.config.clipThresh,
3016 clipFraction=self.config.clipFraction,
3017 )
3018 except RuntimeError as e:
3019 if "Cholesky decomposition failed" in str(e):
3020 raise CholeskyError() from e
3021 else:
3022 raise
3024 self.log.info("WCS fitting done")
3026 outputWcss, cameraParams, colorParams, camera, partialOutputs = self._make_outputs(
3027 wcsf,
3028 inputVisitSummaries,
3029 exposureInfo,
3030 mapTemplate,
3031 extensionInfo,
3032 inputCameraModel=(inputCameraModel if self.config.useInputCameraModel else None),
3033 inputCamera=(inputCamera if self.config.buildCamera else None),
3034 )
3035 outputCatalog = wcsf.getOutputCatalog()
3036 outputCatalog["exposureName"] = np.array(outputCatalog["exposureName"])
3037 outputCatalog["deviceName"] = np.array(outputCatalog["deviceName"])
3039 starCatalog = wcsf.getStarCatalog()
3040 modelParams = self._compute_model_params(wcsf) if self.config.saveModelParams else None
3042 return pipeBase.Struct(
3043 outputWcss=outputWcss,
3044 fitModel=wcsf,
3045 outputCatalog=outputCatalog,
3046 starCatalog=starCatalog,
3047 modelParams=modelParams,
3048 cameraModelParams=cameraParams,
3049 colorParams=colorParams,
3050 camera=camera,
3051 partialOutputs=partialOutputs,
3052 )
3054 def _prep_sky(self, inputVisitSummaries):
3055 """Cluster the input visits into disjoint groups that will define
3056 separate fields in the astrometric fit, and, for each group, get the
3057 convex hull around all of its component visits.
3059 The groups are created such that each visit overlaps with at least one
3060 other visit in the same group by the `visitOverlap` amount, which is
3061 calculated as a fraction of the field-of-view radius, and no visits
3062 from separate groups overlap by more than this amount.
3064 Paramaters
3065 ----------
3066 inputVisitSummaries : `list` [`lsst.afw.table.ExposureCatalog`]
3067 List of catalogs with per-detector summary information.
3069 Returns
3070 -------
3071 fields : `wcsfit.Fields`
3072 Object with field information.
3073 fieldRegions : `dict` [`int`, `lsst.sphgeom.ConvexPolygon`]
3074 Dictionary of regions encompassing each group of input visits,
3075 keyed by an arbitrary index.
3076 """
3077 allDetectorCorners = []
3078 mjds = []
3079 radecs = []
3080 radii = []
3081 for visSum in inputVisitSummaries:
3082 detectorCorners = [
3083 lsst.geom.SpherePoint(ra, dec, lsst.geom.degrees).getVector()
3084 for (ra, dec) in zip(visSum["raCorners"].ravel(), visSum["decCorners"].ravel())
3085 if (np.isfinite(ra) and (np.isfinite(dec)))
3086 ]
3087 if len(detectorCorners) == 0: 3087 ↛ 3091line 3087 didn't jump to line 3091 because the condition on line 3087 was never true
3088 # Skip this visit if none of the detectors have finite ra/dec
3089 # corners, which happens when the WCSs are missing. The visit
3090 # will get formally dropped elsewhere.
3091 continue
3092 allDetectorCorners.append(detectorCorners)
3094 # Get center and approximate radius of field of view
3095 boundingCircle = lsst.sphgeom.ConvexPolygon.convexHull(detectorCorners).getBoundingCircle()
3096 center = lsst.geom.SpherePoint(boundingCircle.getCenter())
3097 ra = center.getRa().asDegrees()
3098 dec = center.getDec().asDegrees()
3099 radecs.append([ra, dec])
3100 radius = boundingCircle.getOpeningAngle()
3101 radii.append(radius)
3103 obsDate = visSum[0].getVisitInfo().getDate()
3104 obsMJD = obsDate.get(obsDate.MJD)
3105 mjds.append(obsMJD)
3107 # Find groups of visits where any one of the visits overlaps another by
3108 # a given fraction of the field-of-view radius.
3109 distance = self.config.visitOverlap * np.median(radii)
3110 clustering = AgglomerativeClustering(
3111 distance_threshold=distance.asDegrees(), n_clusters=None, linkage="single"
3112 )
3113 clusters = clustering.fit(np.array(radecs))
3115 medianMJD = np.median(mjds)
3116 medianEpoch = astropy.time.Time(medianMJD, format="mjd").jyear
3118 fieldNames = []
3119 fieldRAs = []
3120 fieldDecs = []
3121 epochs = []
3122 fieldRegions = {}
3124 for i in range(clusters.n_clusters_):
3125 fieldInd = clusters.labels_ == i
3126 # Concatenate the lists of all detector corners that are in this
3127 # field
3128 fieldDetectors = sum([allDetectorCorners[f] for f, fInd in enumerate(fieldInd) if fInd], [])
3129 hull = lsst.sphgeom.ConvexPolygon.convexHull(fieldDetectors)
3130 center = lsst.geom.SpherePoint(hull.getCentroid())
3131 ra = center.getRa().asDegrees()
3132 dec = center.getDec().asDegrees()
3134 fieldRegions[i] = hull
3135 fieldNames.append(str(i))
3136 fieldRAs.append(ra)
3137 fieldDecs.append(dec)
3138 # Use the same median epoch for all fields so that the final object
3139 # positions are calculated for the same epoch.
3140 epochs.append(medianEpoch)
3142 fields = wcsfit.Fields(fieldNames, fieldRAs, fieldDecs, epochs)
3144 return fields, fieldRegions
3147class GbdesGlobalAstrometricMultibandFitConnections(
3148 GbdesGlobalAstrometricFitConnections,
3149 dimensions=("instrument",),
3150):
3151 outputCatalog = pipeBase.connectionTypes.Output(
3152 doc=(
3153 "Catalog of sources used in fit, along with residuals in pixel coordinates and tangent "
3154 "plane coordinates and chisq values."
3155 ),
3156 name="gbdesGlobalAstrometricMultibandFit_fitStars",
3157 storageClass="ArrowNumpyDict",
3158 dimensions=("instrument",),
3159 )
3160 starCatalog = pipeBase.connectionTypes.Output(
3161 doc=(
3162 "Catalog of best-fit object positions. Also includes the fit proper motion and parallax if "
3163 "fitProperMotion is True."
3164 ),
3165 name="gbdesGlobalAstrometricMultibandFit_starCatalog",
3166 storageClass="ArrowNumpyDict",
3167 dimensions=("instrument",),
3168 )
3169 modelParams = pipeBase.connectionTypes.Output(
3170 doc="WCS parameters and covariance.",
3171 name="gbdesGlobalAstrometricMultibandFit_modelParams",
3172 storageClass="ArrowNumpyDict",
3173 dimensions=("instrument",),
3174 )
3177class GbdesGlobalAstrometricMultibandFitConfig(
3178 GbdesAstrometricFitConfig,
3179 pipelineConnections=GbdesGlobalAstrometricMultibandFitConnections,
3180):
3181 """Configuration for the GbdesGlobalAstrometricMultibandFitTask"""
3183 pass
3186class GbdesGlobalAstrometricMultibandFitTask(GbdesGlobalAstrometricFitTask):
3187 """Calibrate the WCS across multiple visits in multiple filters and
3188 multiple fields using the GBDES package.
3189 """
3191 ConfigClass = GbdesGlobalAstrometricMultibandFitConfig
3192 _DefaultName = "gbdesAstrometricMultibandFit"