LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
jointcal.py
Go to the documentation of this file.
1# This file is part of jointcal.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22import dataclasses
23import collections
24import os
25import logging
26
27import astropy.time
28import numpy as np
29import astropy.units as u
30
31import lsst.geom
32import lsst.utils
33import lsst.pex.config as pexConfig
34import lsst.pipe.base as pipeBase
35from lsst.afw.image import fluxErrFromABMagErr
37import lsst.afw.table
38from lsst.pipe.tasks.colorterms import ColortermLibrary
39from lsst.verify import Job, Measurement
40
41from lsst.meas.algorithms import (ReferenceObjectLoader, ReferenceSourceSelectorTask,
42 LoadReferenceObjectsConfig)
43from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry
44
45import lsst.jointcal
46from lsst.jointcal import MinimizeResult
47
48__all__ = ["JointcalConfig", "JointcalTask"]
49
50Photometry = collections.namedtuple('Photometry', ('fit', 'model'))
51Astrometry = collections.namedtuple('Astrometry', ('fit', 'model', 'sky_to_tan_projection'))
52
53
54# TODO: move this to MeasurementSet in lsst.verify per DM-12655.
55def add_measurement(job, name, value):
56 meas = Measurement(job.metrics[name], value)
57 job.measurements.insert(meas)
58
59
60class JointcalTaskConnections(pipeBase.PipelineTaskConnections,
61 dimensions=("skymap", "tract", "instrument", "physical_filter")):
62 """Middleware input/output connections for jointcal data."""
63 inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
64 doc="The camera instrument that took these observations.",
65 name="camera",
66 storageClass="Camera",
67 dimensions=("instrument",),
68 isCalibration=True,
69 )
70 inputSourceTableVisit = pipeBase.connectionTypes.Input(
71 doc="Source table in parquet format, per visit",
72 name="sourceTable_visit",
73 storageClass="DataFrame",
74 dimensions=("instrument", "visit"),
75 deferLoad=True,
76 multiple=True,
77 )
78 inputVisitSummary = pipeBase.connectionTypes.Input(
79 doc=("Per-visit consolidated exposure metadata built from calexps. "
80 "These catalogs use detector id for the id and must be sorted for "
81 "fast lookups of a detector."),
82 name="visitSummary",
83 storageClass="ExposureCatalog",
84 dimensions=("instrument", "visit"),
85 deferLoad=True,
86 multiple=True,
87 )
88 astrometryRefCat = pipeBase.connectionTypes.PrerequisiteInput(
89 doc="The astrometry reference catalog to match to loaded input catalog sources.",
90 name="gaia_dr2_20200414",
91 storageClass="SimpleCatalog",
92 dimensions=("skypix",),
93 deferLoad=True,
94 multiple=True,
95 )
96 photometryRefCat = pipeBase.connectionTypes.PrerequisiteInput(
97 doc="The photometry reference catalog to match to loaded input catalog sources.",
98 name="ps1_pv3_3pi_20170110",
99 storageClass="SimpleCatalog",
100 dimensions=("skypix",),
101 deferLoad=True,
102 multiple=True,
103 )
104
105 outputWcs = pipeBase.connectionTypes.Output(
106 doc=("Per-tract, per-visit world coordinate systems derived from the fitted model."
107 " These catalogs only contain entries for detectors with an output, and use"
108 " the detector id for the catalog id, sorted on id for fast lookups of a detector."),
109 name="jointcalSkyWcsCatalog",
110 storageClass="ExposureCatalog",
111 dimensions=("instrument", "visit", "skymap", "tract"),
112 multiple=True
113 )
114 outputPhotoCalib = pipeBase.connectionTypes.Output(
115 doc=("Per-tract, per-visit photometric calibrations derived from the fitted model."
116 " These catalogs only contain entries for detectors with an output, and use"
117 " the detector id for the catalog id, sorted on id for fast lookups of a detector."),
118 name="jointcalPhotoCalibCatalog",
119 storageClass="ExposureCatalog",
120 dimensions=("instrument", "visit", "skymap", "tract"),
121 multiple=True
122 )
123
124 # measurements of metrics
125 # The vars() trick used here allows us to set class attributes
126 # programatically. Taken from:
127 # https://stackoverflow.com/questions/2519807/setting-a-class-attribute-with-a-given-name-in-python-while-defining-the-class
128 for name in ("astrometry", "photometry"):
129 vars()[f"{name}_matched_fittedStars"] = pipeBase.connectionTypes.Output(
130 doc=f"The number of cross-matched fittedStars for {name}",
131 name=f"metricvalue_jointcal_{name}_matched_fittedStars",
132 storageClass="MetricValue",
133 dimensions=("skymap", "tract", "instrument", "physical_filter"),
134 )
135 vars()[f"{name}_collected_refStars"] = pipeBase.connectionTypes.Output(
136 doc=f"The number of {name} reference stars drawn from the reference catalog, before matching.",
137 name=f"metricvalue_jointcal_{name}_collected_refStars",
138 storageClass="MetricValue",
139 dimensions=("skymap", "tract", "instrument", "physical_filter"),
140 )
141 vars()[f"{name}_prepared_refStars"] = pipeBase.connectionTypes.Output(
142 doc=f"The number of {name} reference stars matched to fittedStars.",
143 name=f"metricvalue_jointcal_{name}_prepared_refStars",
144 storageClass="MetricValue",
145 dimensions=("skymap", "tract", "instrument", "physical_filter"),
146 )
147 vars()[f"{name}_prepared_fittedStars"] = pipeBase.connectionTypes.Output(
148 doc=f"The number of cross-matched fittedStars after cleanup, for {name}.",
149 name=f"metricvalue_jointcal_{name}_prepared_fittedStars",
150 storageClass="MetricValue",
151 dimensions=("skymap", "tract", "instrument", "physical_filter"),
152 )
153 vars()[f"{name}_prepared_ccdImages"] = pipeBase.connectionTypes.Output(
154 doc=f"The number of ccdImages that will be fit for {name}, after cleanup.",
155 name=f"metricvalue_jointcal_{name}_prepared_ccdImages",
156 storageClass="MetricValue",
157 dimensions=("skymap", "tract", "instrument", "physical_filter"),
158 )
159 vars()[f"{name}_final_chi2"] = pipeBase.connectionTypes.Output(
160 doc=f"The final chi2 of the {name} fit.",
161 name=f"metricvalue_jointcal_{name}_final_chi2",
162 storageClass="MetricValue",
163 dimensions=("skymap", "tract", "instrument", "physical_filter"),
164 )
165 vars()[f"{name}_final_ndof"] = pipeBase.connectionTypes.Output(
166 doc=f"The number of degrees of freedom of the fitted {name}.",
167 name=f"metricvalue_jointcal_{name}_final_ndof",
168 storageClass="MetricValue",
169 dimensions=("skymap", "tract", "instrument", "physical_filter"),
170 )
171
172 def __init__(self, *, config=None):
173 super().__init__(config=config)
174 # When we are only doing one of astrometry or photometry, we don't
175 # need the reference catalog or produce the outputs for the other.
176 # This informs the middleware of that when the QuantumGraph is
177 # generated, so we don't block on getting something we won't need or
178 # create an expectation that downstream tasks will be able to consume
179 # something we won't produce.
180 if not config.doAstrometry:
181 self.prerequisiteInputs.remove("astrometryRefCat")
182 self.outputs.remove("outputWcs")
183 for key in list(self.outputs):
184 if "metricvalue_jointcal_astrometry" in key:
185 self.outputs.remove(key)
186 if not config.doPhotometry:
187 self.prerequisiteInputs.remove("photometryRefCat")
188 self.outputs.remove("outputPhotoCalib")
189 for key in list(self.outputs):
190 if "metricvalue_jointcal_photometry" in key:
191 self.outputs.remove(key)
192
194 return ("inputVisitSummary",)
195
196
197class JointcalConfig(pipeBase.PipelineTaskConfig,
198 pipelineConnections=JointcalTaskConnections):
199 """Configuration for JointcalTask"""
200
201 doAstrometry = pexConfig.Field(
202 doc="Fit astrometry and write the fitted result.",
203 dtype=bool,
204 default=True
205 )
206 doPhotometry = pexConfig.Field(
207 doc="Fit photometry and write the fitted result.",
208 dtype=bool,
209 default=True
210 )
211 sourceFluxType = pexConfig.Field(
212 dtype=str,
213 doc="Source flux field to use in source selection and to get fluxes from the catalog.",
214 default='apFlux_12_0'
215 )
216 positionErrorPedestal = pexConfig.Field(
217 doc="Systematic term to apply to the measured position error (pixels)",
218 dtype=float,
219 default=0.02,
220 )
221 photometryErrorPedestal = pexConfig.Field(
222 doc="Systematic term to apply to the measured error on flux or magnitude as a "
223 "fraction of source flux or magnitude delta (e.g. 0.05 is 5% of flux or +50 millimag).",
224 dtype=float,
225 default=0.0,
226 )
227 # TODO: DM-6885 matchCut should be an geom.Angle
228 matchCut = pexConfig.Field(
229 doc="Matching radius between fitted and reference stars (arcseconds)",
230 dtype=float,
231 default=3.0,
232 )
233 minMeasurements = pexConfig.Field(
234 doc="Minimum number of associated measured stars for a fitted star to be included in the fit",
235 dtype=int,
236 default=2,
237 )
238 minMeasuredStarsPerCcd = pexConfig.Field(
239 doc="Minimum number of measuredStars per ccdImage before printing warnings",
240 dtype=int,
241 default=100,
242 )
243 minRefStarsPerCcd = pexConfig.Field(
244 doc="Minimum number of measuredStars per ccdImage before printing warnings",
245 dtype=int,
246 default=30,
247 )
248 allowLineSearch = pexConfig.Field(
249 doc="Allow a line search during minimization, if it is reasonable for the model"
250 " (models with a significant non-linear component, e.g. constrainedPhotometry).",
251 dtype=bool,
252 default=False
253 )
254 astrometrySimpleOrder = pexConfig.Field(
255 doc="Polynomial order for fitting the simple astrometry model.",
256 dtype=int,
257 default=3,
258 )
259 astrometryChipOrder = pexConfig.Field(
260 doc="Order of the per-chip transform for the constrained astrometry model.",
261 dtype=int,
262 default=1,
263 )
264 astrometryVisitOrder = pexConfig.Field(
265 doc="Order of the per-visit transform for the constrained astrometry model.",
266 dtype=int,
267 default=5,
268 )
269 useInputWcs = pexConfig.Field(
270 doc="Use the input calexp WCSs to initialize a SimpleAstrometryModel.",
271 dtype=bool,
272 default=True,
273 )
274 astrometryModel = pexConfig.ChoiceField(
275 doc="Type of model to fit to astrometry",
276 dtype=str,
277 default="constrained",
278 allowed={"simple": "One polynomial per ccd",
279 "constrained": "One polynomial per ccd, and one polynomial per visit"}
280 )
281 photometryModel = pexConfig.ChoiceField(
282 doc="Type of model to fit to photometry",
283 dtype=str,
284 default="constrainedMagnitude",
285 allowed={"simpleFlux": "One constant zeropoint per ccd and visit, fitting in flux space.",
286 "constrainedFlux": "Constrained zeropoint per ccd, and one polynomial per visit,"
287 " fitting in flux space.",
288 "simpleMagnitude": "One constant zeropoint per ccd and visit,"
289 " fitting in magnitude space.",
290 "constrainedMagnitude": "Constrained zeropoint per ccd, and one polynomial per visit,"
291 " fitting in magnitude space.",
292 }
293 )
294 applyColorTerms = pexConfig.Field(
295 doc="Apply photometric color terms to reference stars?"
296 "Requires that colorterms be set to a ColortermLibrary",
297 dtype=bool,
298 default=False
299 )
300 colorterms = pexConfig.ConfigField(
301 doc="Library of photometric reference catalog name to color term dict.",
302 dtype=ColortermLibrary,
303 )
304 photometryVisitOrder = pexConfig.Field(
305 doc="Order of the per-visit polynomial transform for the constrained photometry model.",
306 dtype=int,
307 default=7,
308 )
309 photometryDoRankUpdate = pexConfig.Field(
310 doc=("Do the rank update step during minimization. "
311 "Skipping this can help deal with models that are too non-linear."),
312 dtype=bool,
313 default=True,
314 )
315 astrometryDoRankUpdate = pexConfig.Field(
316 doc=("Do the rank update step during minimization (should not change the astrometry fit). "
317 "Skipping this can help deal with models that are too non-linear."),
318 dtype=bool,
319 default=True,
320 )
321 outlierRejectSigma = pexConfig.Field(
322 doc="How many sigma to reject outliers at during minimization.",
323 dtype=float,
324 default=5.0,
325 )
326 astrometryOutlierRelativeTolerance = pexConfig.Field(
327 doc=("Convergence tolerance for outlier rejection threshold when fitting astrometry. Iterations will "
328 "stop when the fractional change in the chi2 cut level is below this value. If tolerance is set "
329 "to zero, iterations will continue until there are no more outliers. We suggest a value of 0.002"
330 "as a balance between a shorter minimization runtime and achieving a final fitted model that is"
331 "close to the solution found when removing all outliers."),
332 dtype=float,
333 default=0,
334 )
335 maxPhotometrySteps = pexConfig.Field(
336 doc="Maximum number of minimize iterations to take when fitting photometry.",
337 dtype=int,
338 default=20,
339 )
340 maxAstrometrySteps = pexConfig.Field(
341 doc="Maximum number of minimize iterations to take when fitting astrometry.",
342 dtype=int,
343 default=20,
344 )
345 astrometryRefObjLoader = pexConfig.ConfigField(
346 dtype=LoadReferenceObjectsConfig,
347 doc="Reference object loader for astrometric fit",
348 )
349 photometryRefObjLoader = pexConfig.ConfigField(
350 dtype=LoadReferenceObjectsConfig,
351 doc="Reference object loader for photometric fit",
352 )
353 sourceSelector = sourceSelectorRegistry.makeField(
354 doc="How to select sources for cross-matching",
355 default="science"
356 )
357 astrometryReferenceSelector = pexConfig.ConfigurableField(
358 target=ReferenceSourceSelectorTask,
359 doc="How to down-select the loaded astrometry reference catalog.",
360 )
361 photometryReferenceSelector = pexConfig.ConfigurableField(
362 target=ReferenceSourceSelectorTask,
363 doc="How to down-select the loaded photometry reference catalog.",
364 )
365 astrometryReferenceErr = pexConfig.Field(
366 doc=("Uncertainty on reference catalog coordinates [mas] to use in place of the `coord_*Err` fields. "
367 "If None, then raise an exception if the reference catalog is missing coordinate errors. "
368 "If specified, overrides any existing `coord_*Err` values."),
369 dtype=float,
370 default=None,
371 optional=True
372 )
373
374 # configs for outputting debug information
375 writeInitMatrix = pexConfig.Field(
376 dtype=bool,
377 doc=("Write the pre/post-initialization Hessian and gradient to text files, for debugging. "
378 "Output files will be written to `config.debugOutputPath` and will "
379 "be of the form 'astrometry_[pre|post]init-TRACT-FILTER-mat.txt'. "
380 "Note that these files are the dense versions of the matrix, and so may be very large."),
381 default=False
382 )
383 writeChi2FilesInitialFinal = pexConfig.Field(
384 dtype=bool,
385 doc=("Write .csv files containing the contributions to chi2 for the initialization and final fit. "
386 "Output files will be written to `config.debugOutputPath` and will "
387 "be of the form `astrometry_[initial|final]_chi2-TRACT-FILTER."),
388 default=False
389 )
390 writeChi2FilesOuterLoop = pexConfig.Field(
391 dtype=bool,
392 doc=("Write .csv files containing the contributions to chi2 for the outer fit loop. "
393 "Output files will be written to `config.debugOutputPath` and will "
394 "be of the form `astrometry_init-NN_chi2-TRACT-FILTER`."),
395 default=False
396 )
397 writeInitialModel = pexConfig.Field(
398 dtype=bool,
399 doc=("Write the pre-initialization model to text files, for debugging. "
400 "Output files will be written to `config.debugOutputPath` and will be "
401 "of the form `initial_astrometry_model-TRACT_FILTER.txt`."),
402 default=False
403 )
404 debugOutputPath = pexConfig.Field(
405 dtype=str,
406 default=".",
407 doc=("Path to write debug output files to. Used by "
408 "`writeInitialModel`, `writeChi2Files*`, `writeInitMatrix`.")
409 )
410 detailedProfile = pexConfig.Field(
411 dtype=bool,
412 default=False,
413 doc="Output separate profiling information for different parts of jointcal, e.g. data read, fitting"
414 )
415
416 def validate(self):
417 super().validate()
418 if self.doPhotometry and self.applyColorTerms and len(self.colorterms.data) == 0:
419 msg = "applyColorTerms=True requires the `colorterms` field be set to a ColortermLibrary."
420 raise pexConfig.FieldValidationError(JointcalConfig.colorterms, self, msg)
421 if self.doAstrometry and not self.doPhotometry and self.applyColorTerms:
422 msg = ("Only doing astrometry, but Colorterms are not applied for astrometry;"
423 "applyColorTerms=True will be ignored.")
424 logging.getLogger("lsst.jointcal").warning(msg)
425
426 def setDefaults(self):
427 # Use only primary stars.
428 self.sourceSelector["science"].doRequirePrimary = True
429 # Use only stars because aperture fluxes of galaxies are biased and depend on seeing.
430 self.sourceSelector["science"].doUnresolved = True
431 self.sourceSelector["science"].unresolved.name = "sizeExtendedness"
432 # with dependable signal to noise ratio.
433 self.sourceSelector["science"].doSignalToNoise = True
434 # Min SNR must be > 0 because jointcal cannot handle negative fluxes,
435 # and S/N > 10 to use sources that are not too faint, and thus better measured.
436 self.sourceSelector["science"].signalToNoise.minimum = 10.
437 # Base SNR selection on `sourceFluxType` because that is the flux that jointcal fits.
438 self.sourceSelector["science"].signalToNoise.fluxField = f"{self.sourceFluxType}_instFlux"
439 self.sourceSelector["science"].signalToNoise.errField = f"{self.sourceFluxType}_instFluxErr"
440 # Do not trust blended sources" aperture fluxes which also depend on seeing.
441 self.sourceSelector["science"].doIsolated = True
442 self.sourceSelector["science"].isolated.parentName = "parentSourceId"
443 self.sourceSelector["science"].isolated.nChildName = "deblend_nChild"
444 # Do not trust either flux or centroid measurements with flags,
445 # chosen from the usual QA flags for stars)
446 self.sourceSelector["science"].doFlags = True
447 badFlags = ["pixelFlags_edge",
448 "pixelFlags_saturated",
449 "pixelFlags_interpolatedCenter",
450 "pixelFlags_interpolated",
451 "pixelFlags_crCenter",
452 "pixelFlags_bad",
453 "hsmPsfMoments_flag",
454 f"{self.sourceFluxType}_flag",
455 ]
456 self.sourceSelector["science"].flags.bad = badFlags
457 self.sourceSelector["science"].doRequireFiniteRaDec = True
458 self.sourceSelector["science"].requireFiniteRaDec.raColName = "ra"
459 self.sourceSelector["science"].requireFiniteRaDec.decColName = "dec"
460
461 # Use Gaia-DR2 with proper motions for astrometry; phot_g_mean is the
462 # primary Gaia band, but is not like any normal photometric band.
463 self.astrometryRefObjLoader.requireProperMotion = True
464 self.astrometryRefObjLoader.anyFilterMapsToThis = "phot_g_mean"
465
466
467def writeModel(model, filename, log):
468 """Write model to outfile."""
469 with open(filename, "w") as file:
470 file.write(repr(model))
471 log.info("Wrote %s to file: %s", model, filename)
472
473
474@dataclasses.dataclass
476 """The input data jointcal needs for each detector/visit."""
477 visit: int
478 """The visit identifier of this exposure."""
480 """The catalog derived from this exposure."""
482 """The VisitInfo of this exposure."""
484 """The detector of this exposure."""
486 """The photometric calibration of this exposure."""
488 """The WCS of this exposure."""
490 """The bounding box of this exposure."""
492 """The filter of this exposure."""
493
494
495class JointcalTask(pipeBase.PipelineTask):
496 """Astrometricly and photometricly calibrate across multiple visits of the
497 same field.
498 """
499
500 ConfigClass = JointcalConfig
501 _DefaultName = "jointcal"
502
503 def __init__(self, **kwargs):
504 super().__init__(**kwargs)
505 self.makeSubtask("sourceSelector")
506 if self.config.doAstrometry:
507 self.makeSubtask("astrometryReferenceSelector")
508 else:
510 if self.config.doPhotometry:
511 self.makeSubtask("photometryReferenceSelector")
512 else:
514
515 # To hold various computed metrics for use by tests
516 self.job = Job.load_metrics_package(subset='jointcal')
517
518 def runQuantum(self, butlerQC, inputRefs, outputRefs):
519 # We override runQuantum to set up the refObjLoaders and write the
520 # outputs to the correct refs.
521 inputs = butlerQC.get(inputRefs)
522 # We want the tract number for writing debug files
523 tract = butlerQC.quantum.dataId['tract']
524 if self.config.doAstrometry:
526 dataIds=[ref.datasetRef.dataId
527 for ref in inputRefs.astrometryRefCat],
528 refCats=inputs.pop('astrometryRefCat'),
529 config=self.config.astrometryRefObjLoader,
530 name=self.config.connections.astrometryRefCat,
531 log=self.log)
532 if self.config.doPhotometry:
534 dataIds=[ref.datasetRef.dataId
535 for ref in inputRefs.photometryRefCat],
536 refCats=inputs.pop('photometryRefCat'),
537 config=self.config.photometryRefObjLoader,
538 name=self.config.connections.photometryRefCat,
539 log=self.log)
540 outputs = self.run(**inputs, tract=tract)
541 self._put_metrics(butlerQC, outputs.job, outputRefs)
542 if self.config.doAstrometry:
543 self._put_output(butlerQC, outputs.outputWcs, outputRefs.outputWcs,
544 inputs['inputCamera'], "setWcs")
545 if self.config.doPhotometry:
546 self._put_output(butlerQC, outputs.outputPhotoCalib, outputRefs.outputPhotoCalib,
547 inputs['inputCamera'], "setPhotoCalib")
548
549 def _put_metrics(self, butlerQC, job, outputRefs):
550 """Persist all measured metrics stored in a job.
551
552 Parameters
553 ----------
554 butlerQC : `lsst.pipe.base.QuantumContext`
555 A butler which is specialized to operate in the context of a
556 `lsst.daf.butler.Quantum`; This is the input to `runQuantum`.
557 job : `lsst.verify.job.Job`
558 Measurements of metrics to persist.
559 outputRefs : `list` [`lsst.pipe.base.OutputQuantizedConnection`]
560 The DatasetRefs to persist the data to.
561 """
562 for key in job.measurements.keys():
563 butlerQC.put(job.measurements[key], getattr(outputRefs, key.fqn.replace('jointcal.', '')))
564
565 def _put_output(self, butlerQC, outputs, outputRefs, camera, setter):
566 """Persist the output datasets to their appropriate datarefs.
567
568 Parameters
569 ----------
570 butlerQC : `lsst.pipe.base.QuantumContext`
571 A butler which is specialized to operate in the context of a
572 `lsst.daf.butler.Quantum`; This is the input to `runQuantum`.
573 outputs : `dict` [`tuple`, `lsst.afw.geom.SkyWcs`] or
574 `dict` [`tuple, `lsst.afw.image.PhotoCalib`]
575 The fitted objects to persist.
576 outputRefs : `list` [`lsst.pipe.base.OutputQuantizedConnection`]
577 The DatasetRefs to persist the data to.
578 camera : `lsst.afw.cameraGeom.Camera`
579 The camera for this instrument, to get detector ids from.
580 setter : `str`
581 The method to call on the ExposureCatalog to set each output.
582 """
584 schema.addField('visit', type='L', doc='Visit number')
585
586 def new_catalog(visit, size):
587 """Return an catalog ready to be filled with appropriate output."""
588 catalog = lsst.afw.table.ExposureCatalog(schema)
589 catalog.resize(size)
590 catalog['visit'] = visit
591 metadata = lsst.daf.base.PropertyList()
592 metadata.add("COMMENT", "Catalog id is detector id, sorted.")
593 metadata.add("COMMENT", "Only detectors with data have entries.")
594 return catalog
595
596 # count how many detectors have output for each visit
597 detectors_per_visit = collections.defaultdict(int)
598 for key in outputs:
599 # key is (visit, detector_id), and we only need visit here
600 detectors_per_visit[key[0]] += 1
601
602 for ref in outputRefs:
603 visit = ref.dataId['visit']
604 catalog = new_catalog(visit, detectors_per_visit[visit])
605 # Iterate over every detector and skip the ones we don't have output for.
606 i = 0
607 for detector in camera:
608 detectorId = detector.getId()
609 key = (ref.dataId['visit'], detectorId)
610 if key not in outputs:
611 # skip detectors we don't have output for
612 self.log.debug("No %s output for detector %s in visit %s",
613 setter[3:], detectorId, visit)
614 continue
615
616 catalog[i].setId(detectorId)
617 getattr(catalog[i], setter)(outputs[key])
618 i += 1
619
620 catalog.sort() # ensure that the detectors are in sorted order, for fast lookups
621 butlerQC.put(catalog, ref)
622 self.log.info("Wrote %s detectors to %s", i, ref)
623
624 def run(self, inputSourceTableVisit, inputVisitSummary, inputCamera, tract=None):
625 # Docstring inherited.
626
627 # We take values out of the Parquet table, and put them in "flux_",
628 # and the config.sourceFluxType field is used during that extraction,
629 # so just use "flux" here.
630 sourceFluxField = "flux"
631 jointcalControl = lsst.jointcal.JointcalControl(sourceFluxField)
632 associations = lsst.jointcal.Associations()
633 self.focalPlaneBBox = inputCamera.getFpBBox()
634 oldWcsList, bands = self._load_data(inputSourceTableVisit,
635 inputVisitSummary,
636 associations,
637 jointcalControl,
638 inputCamera)
639
640 boundingCircle, center, radius, defaultFilter, epoch = self._prep_sky(associations, bands)
641
642 if self.config.doAstrometry:
643 astrometry = self._do_load_refcat_and_fit(associations, defaultFilter, center, radius,
644 name="astrometry",
645 refObjLoader=self.astrometryRefObjLoader,
646 referenceSelector=self.astrometryReferenceSelector,
647 fit_function=self._fit_astrometry,
648 tract=tract,
649 epoch=epoch)
650 astrometry_output = self._make_output(associations.getCcdImageList(),
651 astrometry.model,
652 "makeSkyWcs")
653 else:
654 astrometry_output = None
655
656 if self.config.doPhotometry:
657 photometry = self._do_load_refcat_and_fit(associations, defaultFilter, center, radius,
658 name="photometry",
659 refObjLoader=self.photometryRefObjLoader,
660 referenceSelector=self.photometryReferenceSelector,
661 fit_function=self._fit_photometry,
662 tract=tract,
663 epoch=epoch,
664 reject_bad_fluxes=True)
665 photometry_output = self._make_output(associations.getCcdImageList(),
666 photometry.model,
667 "toPhotoCalib")
668 else:
669 photometry_output = None
670
671 return pipeBase.Struct(outputWcs=astrometry_output,
672 outputPhotoCalib=photometry_output,
673 job=self.job,
674 astrometryRefObjLoader=self.astrometryRefObjLoader,
675 photometryRefObjLoader=self.photometryRefObjLoader)
676
677 def _load_data(self, inputSourceTableVisit, inputVisitSummary, associations,
678 jointcalControl, camera):
679 """Read the data that jointcal needs to run.
680
681 Modifies ``associations`` in-place with the loaded data.
682
683 Parameters
684 ----------
685 inputSourceTableVisit : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
686 References to visit-level DataFrames to load the catalog data from.
687 inputVisitSummary : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
688 Visit-level exposure summary catalog with metadata.
689 associations : `lsst.jointcal.Associations`
690 Object to add the loaded data to by constructing new CcdImages.
691 jointcalControl : `jointcal.JointcalControl`
692 Control object for C++ associations management.
693 camera : `lsst.afw.cameraGeom.Camera`
694 Camera object for detector geometry.
695
696 Returns
697 -------
698 oldWcsList: `list` [`lsst.afw.geom.SkyWcs`]
699 The original WCS of the input data, to aid in writing tests.
700 bands : `list` [`str`]
701 The filter bands of each input dataset.
702 """
703 oldWcsList = []
704 filters = []
705 load_cat_profile_file = 'jointcal_load_data.prof' if self.config.detailedProfile else ''
706 with lsst.utils.timer.profile(load_cat_profile_file):
707 table = make_schema_table() # every detector catalog has the same layout
708 # No guarantee that the input is in the same order of visits, so we have to map one of them.
709 catalogMap = {ref.dataId['visit']: i for i, ref in enumerate(inputSourceTableVisit)}
710 detectorDict = {detector.getId(): detector for detector in camera}
711
712 columns = None
713
714 for visitSummaryRef in inputVisitSummary:
715 visitSummary = visitSummaryRef.get()
716
717 dataRef = inputSourceTableVisit[catalogMap[visitSummaryRef.dataId['visit']]]
718 if columns is None:
719 inColumns = dataRef.get(component='columns')
720 columns, ixxColumns = get_sourceTable_visit_columns(inColumns,
721 self.config,
723 visitCatalog = dataRef.get(parameters={'columns': columns})
724
725 selected = self.sourceSelector.run(visitCatalog)
726 if len(selected.sourceCat) == 0:
727 self.log.warning("No sources selected in visit %s. Skipping...",
728 visitSummary["visit"][0])
729 continue
730
731 # Build a CcdImage for each detector in this visit.
732 detectors = {id: index for index, id in enumerate(visitSummary['id'])}
733 for id, index in detectors.items():
735 selected.sourceCat,
736 id,
737 ixxColumns,
738 self.config.sourceFluxType,
739 self.log)
740 if catalog is None:
741 continue
742 data = self._make_one_input_data(visitSummary[index], catalog, detectorDict)
743 result = self._build_ccdImage(data, associations, jointcalControl)
744 if result is not None:
745 oldWcsList.append(result.wcs)
746 # A visit has only one band, so we can just use the first.
747 filters.append(data.filter)
748 if len(filters) == 0:
749 raise RuntimeError("No data to process: did source selector remove all sources?")
750 filters = collections.Counter(filters)
751
752 return oldWcsList, filters
753
754 def _make_one_input_data(self, visitRecord, catalog, detectorDict):
755 """Return a data structure for this detector+visit."""
756 return JointcalInputData(visit=visitRecord['visit'],
757 catalog=catalog,
758 visitInfo=visitRecord.getVisitInfo(),
759 detector=detectorDict[visitRecord.getId()],
760 photoCalib=visitRecord.getPhotoCalib(),
761 wcs=visitRecord.getWcs(),
762 bbox=visitRecord.getBBox(),
763 # ExposureRecord doesn't have a FilterLabel yet,
764 # so we have to make one.
765 filter=lsst.afw.image.FilterLabel(band=visitRecord['band'],
766 physical=visitRecord['physical_filter']))
767
768 def _build_ccdImage(self, data, associations, jointcalControl):
769 """
770 Extract the necessary things from this catalog+metadata to add a new
771 ccdImage.
772
773 Parameters
774 ----------
775 data : `JointcalInputData`
776 The loaded input data.
777 associations : `lsst.jointcal.Associations`
778 Object to add the info to, to construct a new CcdImage
779 jointcalControl : `jointcal.JointcalControl`
780 Control object for associations management
781
782 Returns
783 ------
784 namedtuple or `None`
785 ``wcs``
786 The TAN WCS of this image, read from the calexp
787 (`lsst.afw.geom.SkyWcs`).
788 ``key``
789 A key to identify this dataRef by its visit and ccd ids
790 (`namedtuple`).
791 `None`
792 if there are no sources in the loaded catalog.
793 """
794 if len(data.catalog) == 0:
795 self.log.warning("No sources selected in visit %s ccd %s", data.visit, data.detector.getId())
796 return None
797
798 associations.createCcdImage(data.catalog,
799 data.wcs,
800 data.visitInfo,
801 data.bbox,
802 data.filter.physicalLabel,
803 data.photoCalib,
804 data.detector,
805 data.visit,
806 data.detector.getId(),
807 jointcalControl)
808
809 Result = collections.namedtuple('Result_from_build_CcdImage', ('wcs', 'key'))
810 Key = collections.namedtuple('Key', ('visit', 'ccd'))
811 return Result(data.wcs, Key(data.visit, data.detector.getId()))
812
813 def _getDebugPath(self, filename):
814 """Constructs a path to filename using the configured debug path.
815 """
816 return os.path.join(self.config.debugOutputPath, filename)
817
818 def _prep_sky(self, associations, filters):
819 """Prepare on-sky and other data that must be computed after data has
820 been read.
821 """
822 associations.computeCommonTangentPoint()
823
824 boundingCircle = associations.computeBoundingCircle()
825 center = lsst.geom.SpherePoint(boundingCircle.getCenter())
826 radius = lsst.geom.Angle(boundingCircle.getOpeningAngle().asRadians(), lsst.geom.radians)
827
828 self.log.info(f"Data has center={center} with radius={radius.asDegrees()} degrees.")
829
830 # Determine a default filter band associated with the catalog. See DM-9093
831 defaultFilter = filters.most_common(1)[0][0]
832 self.log.debug("Using '%s' filter for reference flux", defaultFilter.physicalLabel)
833
834 # compute and set the reference epoch of the observations, for proper motion corrections
835 epoch = self._compute_proper_motion_epoch(associations.getCcdImageList())
836 associations.setEpoch(epoch.jyear)
837
838 return boundingCircle, center, radius, defaultFilter, epoch
839
841 """Check whether we should override the refcat coordinate errors, and
842 return the overridden error if necessary.
843
844 Parameters
845 ----------
846 refCat : `lsst.afw.table.SimpleCatalog`
847 The reference catalog to check for a ``coord_raErr`` field.
848 name : `str`
849 Whether we are doing "astrometry" or "photometry".
850
851 Returns
852 -------
853 refCoordErr : `float`
854 The refcat coordinate error to use, or NaN if we are not overriding
855 those fields.
856
857 Raises
858 ------
859 lsst.pex.config.FieldValidationError
860 Raised if the refcat does not contain coordinate errors and
861 ``config.astrometryReferenceErr`` is not set.
862 """
863 # This value doesn't matter for photometry, so just set something to
864 # keep old refcats from causing problems.
865 if name.lower() == "photometry":
866 if 'coord_raErr' not in refCat.schema:
867 return 100
868 else:
869 return float('nan')
870
871 if self.config.astrometryReferenceErr is None and 'coord_raErr' not in refCat.schema:
872 msg = ("Reference catalog does not contain coordinate errors, "
873 "and config.astrometryReferenceErr not supplied.")
874 raise pexConfig.FieldValidationError(JointcalConfig.astrometryReferenceErr,
875 self.config,
876 msg)
877
878 if self.config.astrometryReferenceErr is not None and 'coord_raErr' in refCat.schema:
879 self.log.warning("Overriding reference catalog coordinate errors with %f/coordinate [mas]",
880 self.config.astrometryReferenceErr)
881
882 if self.config.astrometryReferenceErr is None:
883 return float('nan')
884 else:
885 return self.config.astrometryReferenceErr
886
887 def _compute_proper_motion_epoch(self, ccdImageList):
888 """Return the proper motion correction epoch of the provided images.
889
890 Parameters
891 ----------
892 ccdImageList : `list` [`lsst.jointcal.CcdImage`]
893 The images to compute the appropriate epoch for.
894
895 Returns
896 -------
897 epoch : `astropy.time.Time`
898 The date to use for proper motion corrections.
899 """
900 return astropy.time.Time(np.mean([ccdImage.getEpoch() for ccdImage in ccdImageList]),
901 format="jyear",
902 scale="tai")
903
904 def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
905 tract="", match_cut=3.0,
906 reject_bad_fluxes=False, *,
907 name="", refObjLoader=None, referenceSelector=None,
908 fit_function=None, epoch=None):
909 """Load reference catalog, perform the fit, and return the result.
910
911 Parameters
912 ----------
913 associations : `lsst.jointcal.Associations`
914 The star/reference star associations to fit.
915 defaultFilter : `lsst.afw.image.FilterLabel`
916 filter to load from reference catalog.
917 center : `lsst.geom.SpherePoint`
918 ICRS center of field to load from reference catalog.
919 radius : `lsst.geom.Angle`
920 On-sky radius to load from reference catalog.
921 name : `str`
922 Name of thing being fit: "astrometry" or "photometry".
923 refObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`
924 Reference object loader to use to load a reference catalog.
925 referenceSelector : `lsst.meas.algorithms.ReferenceSourceSelectorTask`
926 Selector to use to pick objects from the loaded reference catalog.
927 fit_function : callable
928 Function to call to perform fit (takes Associations object).
929 tract : `str`, optional
930 Name of tract currently being fit.
931 match_cut : `float`, optional
932 Radius in arcseconds to find cross-catalog matches to during
933 associations.associateCatalogs.
934 reject_bad_fluxes : `bool`, optional
935 Reject refCat sources with NaN/inf flux or NaN/0 fluxErr.
936 epoch : `astropy.time.Time`, optional
937 Epoch to which to correct refcat proper motion and parallax,
938 or `None` to not apply such corrections.
939
940 Returns
941 -------
942 result : `Photometry` or `Astrometry`
943 Result of `fit_function()`
944 """
945 self.log.info("====== Now processing %s...", name)
946 # TODO: this should not print "trying to invert a singular transformation:"
947 # if it does that, something's not right about the WCS...
948 associations.associateCatalogs(match_cut)
949 add_measurement(self.job, 'jointcal.%s_matched_fittedStars' % name,
950 associations.fittedStarListSize())
951
952 applyColorterms = False if name.lower() == "astrometry" else self.config.applyColorTerms
953 refCat, fluxField = self._load_reference_catalog(refObjLoader, referenceSelector,
954 center, radius, defaultFilter,
955 applyColorterms=applyColorterms,
956 epoch=epoch)
957 refCoordErr = self._get_refcat_coordinate_error_override(refCat, name)
958
959 associations.collectRefStars(refCat,
960 self.config.matchCut*lsst.geom.arcseconds,
961 fluxField,
962 refCoordinateErr=refCoordErr,
963 rejectBadFluxes=reject_bad_fluxes)
964 add_measurement(self.job, 'jointcal.%s_collected_refStars' % name,
965 associations.refStarListSize())
966
967 associations.prepareFittedStars(self.config.minMeasurements)
968
969 self._check_star_lists(associations, name)
970 add_measurement(self.job, 'jointcal.%s_prepared_refStars' % name,
971 associations.nFittedStarsWithAssociatedRefStar())
972 add_measurement(self.job, 'jointcal.%s_prepared_fittedStars' % name,
973 associations.fittedStarListSize())
974 add_measurement(self.job, 'jointcal.%s_prepared_ccdImages' % name,
975 associations.nCcdImagesValidForFit())
976
977 fit_profile_file = 'jointcal_fit_%s.prof'%name if self.config.detailedProfile else ''
978 dataName = "{}_{}".format(tract, defaultFilter.physicalLabel)
979 with lsst.utils.timer.profile(fit_profile_file):
980 result = fit_function(associations, dataName)
981 # TODO DM-12446: turn this into a "butler save" somehow.
982 # Save reference and measurement chi2 contributions for this data
983 if self.config.writeChi2FilesInitialFinal:
984 baseName = self._getDebugPath(f"{name}_final_chi2-{dataName}")
985 result.fit.saveChi2Contributions(baseName+"{type}")
986 self.log.info("Wrote chi2 contributions files: %s", baseName)
987
988 return result
989
990 def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radius, filterLabel,
991 applyColorterms=False, epoch=None):
992 """Load the necessary reference catalog sources, convert fluxes to
993 correct units, and apply color term corrections if requested.
994
995 Parameters
996 ----------
997 refObjLoader : `lsst.meas.algorithms.ReferenceObjectLoader`
998 The reference catalog loader to use to get the data.
999 referenceSelector : `lsst.meas.algorithms.ReferenceSourceSelectorTask`
1000 Source selector to apply to loaded reference catalog.
1001 center : `lsst.geom.SpherePoint`
1002 The center around which to load sources.
1003 radius : `lsst.geom.Angle`
1004 The radius around ``center`` to load sources in.
1005 filterLabel : `lsst.afw.image.FilterLabel`
1006 The camera filter to load fluxes for.
1007 applyColorterms : `bool`
1008 Apply colorterm corrections to the refcat for ``filterName``?
1009 epoch : `astropy.time.Time`, optional
1010 Epoch to which to correct refcat proper motion and parallax,
1011 or `None` to not apply such corrections.
1012
1013 Returns
1014 -------
1015 refCat : `lsst.afw.table.SimpleCatalog`
1016 The loaded reference catalog.
1017 fluxField : `str`
1018 The name of the reference catalog flux field appropriate for ``filterName``.
1019 """
1020 skyCircle = refObjLoader.loadSkyCircle(center,
1021 radius,
1022 filterLabel.bandLabel,
1023 epoch=epoch)
1024
1025 selected = referenceSelector.run(skyCircle.refCat)
1026 # Need memory contiguity to get reference filters as a vector.
1027 if not selected.sourceCat.isContiguous():
1028 refCat = selected.sourceCat.copy(deep=True)
1029 else:
1030 refCat = selected.sourceCat
1031
1032 if applyColorterms:
1033 refCatName = refObjLoader.name
1034 self.log.info("Applying color terms for physical filter=%r reference catalog=%s",
1035 filterLabel.physicalLabel, refCatName)
1036 colorterm = self.config.colorterms.getColorterm(filterLabel.physicalLabel,
1037 refCatName,
1038 doRaise=True)
1039
1040 refMag, refMagErr = colorterm.getCorrectedMagnitudes(refCat)
1041 refCat[skyCircle.fluxField] = u.Magnitude(refMag, u.ABmag).to_value(u.nJy)
1042 # TODO: I didn't want to use this, but I'll deal with it in DM-16903
1043 refCat[skyCircle.fluxField+'Err'] = fluxErrFromABMagErr(refMagErr, refMag) * 1e9
1044
1045 return refCat, skyCircle.fluxField
1046
1047 def _check_star_lists(self, associations, name):
1048 # TODO: these should be len(blah), but we need this properly wrapped first.
1049 if associations.nCcdImagesValidForFit() == 0:
1050 raise RuntimeError('No images in the ccdImageList!')
1051 if associations.fittedStarListSize() == 0:
1052 raise RuntimeError('No stars in the {} fittedStarList!'.format(name))
1053 if associations.refStarListSize() == 0:
1054 raise RuntimeError('No stars in the {} reference star list!'.format(name))
1055
1056 def _logChi2AndValidate(self, associations, fit, model, chi2Label, writeChi2Name=None):
1057 """Compute chi2, log it, validate the model, and return chi2.
1058
1059 Parameters
1060 ----------
1061 associations : `lsst.jointcal.Associations`
1062 The star/reference star associations to fit.
1063 fit : `lsst.jointcal.FitterBase`
1064 The fitter to use for minimization.
1065 model : `lsst.jointcal.Model`
1066 The model being fit.
1067 chi2Label : `str`
1068 Label to describe the chi2 (e.g. "Initialized", "Final").
1069 writeChi2Name : `str`, optional
1070 Filename prefix to write the chi2 contributions to.
1071 Do not supply an extension: an appropriate one will be added.
1072
1073 Returns
1074 -------
1075 chi2: `lsst.jointcal.Chi2Accumulator`
1076 The chi2 object for the current fitter and model.
1077
1078 Raises
1079 ------
1080 FloatingPointError
1081 Raised if chi2 is infinite or NaN.
1082 ValueError
1083 Raised if the model is not valid.
1084 """
1085 if writeChi2Name is not None:
1086 fullpath = self._getDebugPath(writeChi2Name)
1087 fit.saveChi2Contributions(fullpath+"{type}")
1088 self.log.info("Wrote chi2 contributions files: %s", fullpath)
1089
1090 chi2 = fit.computeChi2()
1091 self.log.info("%s %s", chi2Label, chi2)
1092 self._check_stars(associations)
1093 if not np.isfinite(chi2.chi2):
1094 raise FloatingPointError(f'{chi2Label} chi2 is invalid: {chi2}')
1095 if not model.validate(associations.getCcdImageList(), chi2.ndof):
1096 raise ValueError("Model is not valid: check log messages for warnings.")
1097 return chi2
1098
1099 def _fit_photometry(self, associations, dataName=None):
1100 """
1101 Fit the photometric data.
1102
1103 Parameters
1104 ----------
1105 associations : `lsst.jointcal.Associations`
1106 The star/reference star associations to fit.
1107 dataName : `str`
1108 Name of the data being processed (e.g. "1234_HSC-Y"), for
1109 identifying debugging files.
1110
1111 Returns
1112 -------
1113 fit_result : `namedtuple`
1114 fit : `lsst.jointcal.PhotometryFit`
1115 The photometric fitter used to perform the fit.
1116 model : `lsst.jointcal.PhotometryModel`
1117 The photometric model that was fit.
1118 """
1119 self.log.info("=== Starting photometric fitting...")
1120
1121 # TODO: should use pex.config.RegistryField here (see DM-9195)
1122 if self.config.photometryModel == "constrainedFlux":
1123 model = lsst.jointcal.ConstrainedFluxModel(associations.getCcdImageList(),
1124 self.focalPlaneBBox,
1125 visitOrder=self.config.photometryVisitOrder,
1126 errorPedestal=self.config.photometryErrorPedestal)
1127 # potentially nonlinear problem, so we may need a line search to converge.
1128 doLineSearch = self.config.allowLineSearch
1129 elif self.config.photometryModel == "constrainedMagnitude":
1130 model = lsst.jointcal.ConstrainedMagnitudeModel(associations.getCcdImageList(),
1131 self.focalPlaneBBox,
1132 visitOrder=self.config.photometryVisitOrder,
1133 errorPedestal=self.config.photometryErrorPedestal)
1134 # potentially nonlinear problem, so we may need a line search to converge.
1135 doLineSearch = self.config.allowLineSearch
1136 elif self.config.photometryModel == "simpleFlux":
1137 model = lsst.jointcal.SimpleFluxModel(associations.getCcdImageList(),
1138 errorPedestal=self.config.photometryErrorPedestal)
1139 doLineSearch = False # purely linear in model parameters, so no line search needed
1140 elif self.config.photometryModel == "simpleMagnitude":
1141 model = lsst.jointcal.SimpleMagnitudeModel(associations.getCcdImageList(),
1142 errorPedestal=self.config.photometryErrorPedestal)
1143 doLineSearch = False # purely linear in model parameters, so no line search needed
1144
1145 fit = lsst.jointcal.PhotometryFit(associations, model)
1146 # TODO DM-12446: turn this into a "butler save" somehow.
1147 # Save reference and measurement chi2 contributions for this data
1148 if self.config.writeChi2FilesInitialFinal:
1149 baseName = f"photometry_initial_chi2-{dataName}"
1150 else:
1151 baseName = None
1152 if self.config.writeInitialModel:
1153 fullpath = self._getDebugPath(f"initial_photometry_model-{dataName}.txt")
1154 writeModel(model, fullpath, self.log)
1155 self._logChi2AndValidate(associations, fit, model, "Initialized", writeChi2Name=baseName)
1156
1157 def getChi2Name(whatToFit):
1158 if self.config.writeChi2FilesOuterLoop:
1159 return f"photometry_init-%s_chi2-{dataName}" % whatToFit
1160 else:
1161 return None
1162
1163 # The constrained model needs the visit transform fit first; the chip
1164 # transform is initialized from the singleFrame PhotoCalib, so it's close.
1165 if self.config.writeInitMatrix:
1166 dumpMatrixFile = self._getDebugPath(f"photometry_preinit-{dataName}")
1167 else:
1168 dumpMatrixFile = ""
1169 if self.config.photometryModel.startswith("constrained"):
1170 # no line search: should be purely (or nearly) linear,
1171 # and we want a large step size to initialize with.
1172 fit.minimize("ModelVisit", dumpMatrixFile=dumpMatrixFile)
1173 self._logChi2AndValidate(associations, fit, model, "Initialize ModelVisit",
1174 writeChi2Name=getChi2Name("ModelVisit"))
1175 dumpMatrixFile = "" # so we don't redo the output on the next step
1176
1177 fit.minimize("Model", doLineSearch=doLineSearch, dumpMatrixFile=dumpMatrixFile)
1178 self._logChi2AndValidate(associations, fit, model, "Initialize Model",
1179 writeChi2Name=getChi2Name("Model"))
1180
1181 fit.minimize("Fluxes") # no line search: always purely linear.
1182 self._logChi2AndValidate(associations, fit, model, "Initialize Fluxes",
1183 writeChi2Name=getChi2Name("Fluxes"))
1184
1185 fit.minimize("Model Fluxes", doLineSearch=doLineSearch)
1186 self._logChi2AndValidate(associations, fit, model, "Initialize ModelFluxes",
1187 writeChi2Name=getChi2Name("ModelFluxes"))
1188
1189 model.freezeErrorTransform()
1190 self.log.debug("Photometry error scales are frozen.")
1191
1192 chi2 = self._iterate_fit(associations,
1193 fit,
1194 self.config.maxPhotometrySteps,
1195 "photometry",
1196 "Model Fluxes",
1197 doRankUpdate=self.config.photometryDoRankUpdate,
1198 doLineSearch=doLineSearch,
1199 dataName=dataName)
1200
1201 add_measurement(self.job, 'jointcal.photometry_final_chi2', chi2.chi2)
1202 add_measurement(self.job, 'jointcal.photometry_final_ndof', chi2.ndof)
1203 return Photometry(fit, model)
1204
1205 def _fit_astrometry(self, associations, dataName=None):
1206 """
1207 Fit the astrometric data.
1208
1209 Parameters
1210 ----------
1211 associations : `lsst.jointcal.Associations`
1212 The star/reference star associations to fit.
1213 dataName : `str`
1214 Name of the data being processed (e.g. "1234_HSC-Y"), for
1215 identifying debugging files.
1216
1217 Returns
1218 -------
1219 fit_result : `namedtuple`
1220 fit : `lsst.jointcal.AstrometryFit`
1221 The astrometric fitter used to perform the fit.
1222 model : `lsst.jointcal.AstrometryModel`
1223 The astrometric model that was fit.
1224 sky_to_tan_projection : `lsst.jointcal.ProjectionHandler`
1225 The model for the sky to tangent plane projection that was used in the fit.
1226 """
1227
1228 self.log.info("=== Starting astrometric fitting...")
1229
1230 associations.deprojectFittedStars()
1231
1232 # NOTE: need to return sky_to_tan_projection so that it doesn't get garbage collected.
1233 # TODO: could we package sky_to_tan_projection and model together so we don't have to manage
1234 # them so carefully?
1235 sky_to_tan_projection = lsst.jointcal.OneTPPerVisitHandler(associations.getCcdImageList())
1236
1237 if self.config.astrometryModel == "constrained":
1238 model = lsst.jointcal.ConstrainedAstrometryModel(associations.getCcdImageList(),
1239 sky_to_tan_projection,
1240 chipOrder=self.config.astrometryChipOrder,
1241 visitOrder=self.config.astrometryVisitOrder)
1242 elif self.config.astrometryModel == "simple":
1243 model = lsst.jointcal.SimpleAstrometryModel(associations.getCcdImageList(),
1244 sky_to_tan_projection,
1245 self.config.useInputWcs,
1246 nNotFit=0,
1247 order=self.config.astrometrySimpleOrder)
1248
1249 fit = lsst.jointcal.AstrometryFit(associations, model, self.config.positionErrorPedestal)
1250 # TODO DM-12446: turn this into a "butler save" somehow.
1251 # Save reference and measurement chi2 contributions for this data
1252 if self.config.writeChi2FilesInitialFinal:
1253 baseName = f"astrometry_initial_chi2-{dataName}"
1254 else:
1255 baseName = None
1256 if self.config.writeInitialModel:
1257 fullpath = self._getDebugPath(f"initial_astrometry_model-{dataName}.txt")
1258 writeModel(model, fullpath, self.log)
1259 self._logChi2AndValidate(associations, fit, model, "Initial", writeChi2Name=baseName)
1260
1261 def getChi2Name(whatToFit):
1262 if self.config.writeChi2FilesOuterLoop:
1263 return f"astrometry_init-%s_chi2-{dataName}" % whatToFit
1264 else:
1265 return None
1266
1267 if self.config.writeInitMatrix:
1268 dumpMatrixFile = self._getDebugPath(f"astrometry_preinit-{dataName}")
1269 else:
1270 dumpMatrixFile = ""
1271 # The constrained model needs the visit transform fit first; the chip
1272 # transform is initialized from the detector's cameraGeom, so it's close.
1273 if self.config.astrometryModel == "constrained":
1274 fit.minimize("DistortionsVisit", dumpMatrixFile=dumpMatrixFile)
1275 self._logChi2AndValidate(associations, fit, model, "Initialize DistortionsVisit",
1276 writeChi2Name=getChi2Name("DistortionsVisit"))
1277 dumpMatrixFile = "" # so we don't redo the output on the next step
1278
1279 fit.minimize("Distortions", dumpMatrixFile=dumpMatrixFile)
1280 self._logChi2AndValidate(associations, fit, model, "Initialize Distortions",
1281 writeChi2Name=getChi2Name("Distortions"))
1282
1283 fit.minimize("Positions")
1284 self._logChi2AndValidate(associations, fit, model, "Initialize Positions",
1285 writeChi2Name=getChi2Name("Positions"))
1286
1287 fit.minimize("Distortions Positions")
1288 self._logChi2AndValidate(associations, fit, model, "Initialize DistortionsPositions",
1289 writeChi2Name=getChi2Name("DistortionsPositions"))
1290
1291 chi2 = self._iterate_fit(associations,
1292 fit,
1293 self.config.maxAstrometrySteps,
1294 "astrometry",
1295 "Distortions Positions",
1296 sigmaRelativeTolerance=self.config.astrometryOutlierRelativeTolerance,
1297 doRankUpdate=self.config.astrometryDoRankUpdate,
1298 dataName=dataName)
1299
1300 add_measurement(self.job, 'jointcal.astrometry_final_chi2', chi2.chi2)
1301 add_measurement(self.job, 'jointcal.astrometry_final_ndof', chi2.ndof)
1302
1303 return Astrometry(fit, model, sky_to_tan_projection)
1304
1305 def _check_stars(self, associations):
1306 """Count measured and reference stars per ccd and warn/log them."""
1307 for ccdImage in associations.getCcdImageList():
1308 nMeasuredStars, nRefStars = ccdImage.countStars()
1309 self.log.debug("ccdImage %s has %s measured and %s reference stars",
1310 ccdImage.getName(), nMeasuredStars, nRefStars)
1311 if nMeasuredStars < self.config.minMeasuredStarsPerCcd:
1312 self.log.warning("ccdImage %s has only %s measuredStars (desired %s)",
1313 ccdImage.getName(), nMeasuredStars, self.config.minMeasuredStarsPerCcd)
1314 if nRefStars < self.config.minRefStarsPerCcd:
1315 self.log.warning("ccdImage %s has only %s RefStars (desired %s)",
1316 ccdImage.getName(), nRefStars, self.config.minRefStarsPerCcd)
1317
1318 def _iterate_fit(self, associations, fitter, max_steps, name, whatToFit,
1319 dataName="",
1320 sigmaRelativeTolerance=0,
1321 doRankUpdate=True,
1322 doLineSearch=False):
1323 """Run fitter.minimize up to max_steps times, returning the final chi2.
1324
1325 Parameters
1326 ----------
1327 associations : `lsst.jointcal.Associations`
1328 The star/reference star associations to fit.
1329 fitter : `lsst.jointcal.FitterBase`
1330 The fitter to use for minimization.
1331 max_steps : `int`
1332 Maximum number of steps to run outlier rejection before declaring
1333 convergence failure.
1334 name : {'photometry' or 'astrometry'}
1335 What type of data are we fitting (for logs and debugging files).
1336 whatToFit : `str`
1337 Passed to ``fitter.minimize()`` to define the parameters to fit.
1338 dataName : `str`, optional
1339 Descriptive name for this dataset (e.g. tract and filter),
1340 for debugging.
1341 sigmaRelativeTolerance : `float`, optional
1342 Convergence tolerance for the fractional change in the chi2 cut
1343 level for determining outliers. If set to zero, iterations will
1344 continue until there are no outliers.
1345 doRankUpdate : `bool`, optional
1346 Do an Eigen rank update during minimization, or recompute the full
1347 matrix and gradient?
1348 doLineSearch : `bool`, optional
1349 Do a line search for the optimum step during minimization?
1350
1351 Returns
1352 -------
1353 chi2: `lsst.jointcal.Chi2Statistic`
1354 The final chi2 after the fit converges, or is forced to end.
1355
1356 Raises
1357 ------
1358 FloatingPointError
1359 Raised if the fitter fails with a non-finite value.
1360 RuntimeError
1361 Raised if the fitter fails for some other reason;
1362 log messages will provide further details.
1363 """
1364 if self.config.writeInitMatrix:
1365 dumpMatrixFile = self._getDebugPath(f"{name}_postinit-{dataName}")
1366 else:
1367 dumpMatrixFile = ""
1368 oldChi2 = lsst.jointcal.Chi2Statistic()
1369 oldChi2.chi2 = float("inf")
1370 for i in range(max_steps):
1371 if self.config.writeChi2FilesOuterLoop:
1372 writeChi2Name = f"{name}_iterate_{i}_chi2-{dataName}"
1373 else:
1374 writeChi2Name = None
1375 result = fitter.minimize(whatToFit,
1376 self.config.outlierRejectSigma,
1377 sigmaRelativeTolerance=sigmaRelativeTolerance,
1378 doRankUpdate=doRankUpdate,
1379 doLineSearch=doLineSearch,
1380 dumpMatrixFile=dumpMatrixFile)
1381 dumpMatrixFile = "" # clear it so we don't write the matrix again.
1382 chi2 = self._logChi2AndValidate(associations, fitter, fitter.getModel(),
1383 f"Fit iteration {i}", writeChi2Name=writeChi2Name)
1384
1385 if result == MinimizeResult.Converged:
1386 if doRankUpdate:
1387 self.log.debug("fit has converged - no more outliers - redo minimization "
1388 "one more time in case we have lost accuracy in rank update.")
1389 # Redo minimization one more time in case we have lost accuracy in rank update
1390 result = fitter.minimize(whatToFit, self.config.outlierRejectSigma,
1391 sigmaRelativeTolerance=sigmaRelativeTolerance)
1392 chi2 = self._logChi2AndValidate(associations, fitter, fitter.getModel(), "Fit completed")
1393
1394 # log a message for a large final chi2, TODO: DM-15247 for something better
1395 if chi2.chi2/chi2.ndof >= 4.0:
1396 self.log.error("Potentially bad fit: High chi-squared/ndof.")
1397
1398 break
1399 elif result == MinimizeResult.Chi2Increased:
1400 self.log.warning("Still some outliers remaining but chi2 increased - retry")
1401 # Check whether the increase was large enough to cause trouble.
1402 chi2Ratio = chi2.chi2 / oldChi2.chi2
1403 if chi2Ratio > 1.5:
1404 self.log.warning('Significant chi2 increase by a factor of %.4g / %.4g = %.4g',
1405 chi2.chi2, oldChi2.chi2, chi2Ratio)
1406 # Based on a variety of HSC jointcal logs (see DM-25779), it
1407 # appears that chi2 increases more than a factor of ~2 always
1408 # result in the fit diverging rapidly and ending at chi2 > 1e10.
1409 # Using 10 as the "failure" threshold gives some room between
1410 # leaving a warning and bailing early.
1411 if chi2Ratio > 10:
1412 msg = ("Large chi2 increase between steps: fit likely cannot converge."
1413 " Try setting one or more of the `writeChi2*` config fields and looking"
1414 " at how individual star chi2-values evolve during the fit.")
1415 raise RuntimeError(msg)
1416 oldChi2 = chi2
1417 elif result == MinimizeResult.NonFinite:
1418 filename = self._getDebugPath("{}_failure-nonfinite_chi2-{}.csv".format(name, dataName))
1419 # TODO DM-12446: turn this into a "butler save" somehow.
1420 fitter.saveChi2Contributions(filename+"{type}")
1421 msg = "Nonfinite value in chi2 minimization, cannot complete fit. Dumped star tables to: {}"
1422 raise FloatingPointError(msg.format(filename))
1423 elif result == MinimizeResult.Failed:
1424 raise RuntimeError("Chi2 minimization failure, cannot complete fit.")
1425 else:
1426 raise RuntimeError("Unxepected return code from minimize().")
1427 else:
1428 self.log.error("%s failed to converge after %d steps"%(name, max_steps))
1429
1430 return chi2
1431
1432 def _make_output(self, ccdImageList, model, func):
1433 """Return the internal jointcal models converted to the afw
1434 structures that will be saved to disk.
1435
1436 Parameters
1437 ----------
1438 ccdImageList : `lsst.jointcal.CcdImageList`
1439 The list of CcdImages to get the output for.
1440 model : `lsst.jointcal.AstrometryModel` or `lsst.jointcal.PhotometryModel`
1441 The internal jointcal model to convert for each `lsst.jointcal.CcdImage`.
1442 func : `str`
1443 The name of the function to call on ``model`` to get the converted
1444 structure. Must accept an `lsst.jointcal.CcdImage`.
1445
1446 Returns
1447 -------
1448 output : `dict` [`tuple`, `lsst.jointcal.AstrometryModel`] or
1449 `dict` [`tuple`, `lsst.jointcal.PhotometryModel`]
1450 The data to be saved, keyed on (visit, detector).
1451 """
1452 output = {}
1453 for ccdImage in ccdImageList:
1454 ccd = ccdImage.ccdId
1455 visit = ccdImage.visit
1456 self.log.debug("%s for visit: %d, ccd: %d", func, visit, ccd)
1457 output[(visit, ccd)] = getattr(model, func)(ccdImage)
1458 return output
1459
1460
1462 """Return an afw SourceTable to use as a base for creating the
1463 SourceCatalog to insert values from the dataFrame into.
1464
1465 Returns
1466 -------
1467 table : `lsst.afw.table.SourceTable`
1468 Table with schema and slots to use to make SourceCatalogs.
1469 """
1471 schema.addField("centroid_x", "D")
1472 schema.addField("centroid_y", "D")
1473 schema.addField("centroid_xErr", "F")
1474 schema.addField("centroid_yErr", "F")
1475 schema.addField("shape_xx", "D")
1476 schema.addField("shape_yy", "D")
1477 schema.addField("shape_xy", "D")
1478 schema.addField("flux_instFlux", "D")
1479 schema.addField("flux_instFluxErr", "D")
1480 table = lsst.afw.table.SourceTable.make(schema)
1481 table.defineCentroid("centroid")
1482 table.defineShape("shape")
1483 return table
1484
1485
1486def get_sourceTable_visit_columns(inColumns, config, sourceSelector):
1487 """
1488 Get the sourceTable_visit columns to load from the catalogs.
1489
1490 Parameters
1491 ----------
1492 inColumns : `list`
1493 List of columns known to be available in the sourceTable_visit.
1494 config : `JointcalConfig`
1495 A filled-in config to to help define column names.
1496 sourceSelector : `lsst.meas.algorithms.BaseSourceSelectorTask`
1497 A configured source selector to define column names to load.
1498
1499 Returns
1500 -------
1501 columns : `list`
1502 List of columns to read from sourceTable_visit.
1503 ixxColumns : `list`
1504 Name of the ixx/iyy/ixy columns.
1505 """
1506 columns = ['visit', 'detector',
1507 'sourceId', 'x', 'xErr', 'y', 'yErr',
1508 config.sourceFluxType + '_instFlux', config.sourceFluxType + '_instFluxErr']
1509
1510 if 'ixx' in inColumns:
1511 # New columns post-DM-31825
1512 ixxColumns = ['ixx', 'iyy', 'ixy']
1513 else:
1514 # Old columns pre-DM-31825
1515 ixxColumns = ['Ixx', 'Iyy', 'Ixy']
1516 columns.extend(ixxColumns)
1517
1518 if sourceSelector.config.doFlags:
1519 columns.extend(sourceSelector.config.flags.bad)
1520 if sourceSelector.config.doUnresolved:
1521 columns.append(sourceSelector.config.unresolved.name)
1522 if sourceSelector.config.doIsolated:
1523 columns.append(sourceSelector.config.isolated.parentName)
1524 columns.append(sourceSelector.config.isolated.nChildName)
1525 if sourceSelector.config.doRequireFiniteRaDec:
1526 columns.append(sourceSelector.config.requireFiniteRaDec.raColName)
1527 columns.append(sourceSelector.config.requireFiniteRaDec.decColName)
1528 if sourceSelector.config.doRequirePrimary:
1529 columns.append(sourceSelector.config.requirePrimary.primaryColName)
1530
1531 return columns, ixxColumns
1532
1533
1534def extract_detector_catalog_from_visit_catalog(table, visitCatalog, detectorId,
1535 ixxColumns, sourceFluxType, log):
1536 """Return an afw SourceCatalog extracted from a visit-level dataframe,
1537 limited to just one detector.
1538
1539 Parameters
1540 ----------
1541 table : `lsst.afw.table.SourceTable`
1542 Table factory to use to make the SourceCatalog that will be
1543 populated with data from ``visitCatalog``.
1544 visitCatalog : `pandas.DataFrame`
1545 DataFrame to extract a detector catalog from.
1546 detectorId : `int`
1547 Numeric id of the detector to extract from ``visitCatalog``.
1548 ixxColumns : `list` [`str`]
1549 Names of the ixx/iyy/ixy columns in the catalog.
1550 sourceFluxType : `str`
1551 Name of the catalog field to load instFluxes from.
1552 log : `logging.Logger`
1553 Logging instance to log to.
1554
1555 Returns
1556 -------
1557 catalog : `lsst.afw.table.SourceCatalog`, or `None`
1558 Detector-level catalog extracted from ``visitCatalog``, or `None`
1559 if there was no data to load.
1560 """
1561 # map from dataFrame column to afw table column
1562 mapping = {'x': 'centroid_x',
1563 'y': 'centroid_y',
1564 'xErr': 'centroid_xErr',
1565 'yErr': 'centroid_yErr',
1566 ixxColumns[0]: 'shape_xx',
1567 ixxColumns[1]: 'shape_yy',
1568 ixxColumns[2]: 'shape_xy',
1569 f'{sourceFluxType}_instFlux': 'flux_instFlux',
1570 f'{sourceFluxType}_instFluxErr': 'flux_instFluxErr',
1571 }
1572
1573 catalog = lsst.afw.table.SourceCatalog(table)
1574 matched = visitCatalog['detector'] == detectorId
1575 n = sum(matched)
1576 if n == 0:
1577 return None
1578 catalog.resize(sum(matched))
1579 view = visitCatalog.loc[matched]
1580 catalog['id'] = view.index
1581 for dfCol, afwCol in mapping.items():
1582 catalog[afwCol] = view[dfCol]
1583
1584 log.debug("%d sources selected in visit %d detector %d",
1585 len(catalog),
1586 view['visit'].iloc[0], # all visits in this catalog are the same, so take the first
1587 detectorId)
1588 return catalog
A representation of a detector in a mosaic camera.
Definition Detector.h:185
A group of labels for a filter in an exposure or coadd.
Definition FilterLabel.h:58
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
Information about a single exposure of an imaging camera.
Definition VisitInfo.h:68
Custom catalog class for ExposureRecord/Table.
Definition Exposure.h:311
static Schema makeMinimalSchema()
Return a minimal schema for Exposure tables and records.
Definition Exposure.h:216
static std::shared_ptr< SourceTable > make(Schema const &schema, std::shared_ptr< IdFactory > const &idFactory)
Construct a new table.
Definition Source.cc:400
static Schema makeMinimalSchema()
Return a minimal schema for Source tables and records.
Definition Source.h:258
Class for storing ordered metadata with comments.
A class representing an angle.
Definition Angle.h:128
An integer coordinate rectangle.
Definition Box.h:55
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
The class that implements the relations between MeasuredStar and FittedStar.
Class that handles the astrometric least squares problem.
Simple structure to accumulate chi2 and ndof.
Definition Chi2.h:52
A multi-component model, fitting mappings for sensors and visits simultaneously.
A projection handler in which all CCDs from the same visit have the same tangent point.
Class that handles the photometric least squares problem.
A model where there is one independent transform per CcdImage.
_build_ccdImage(self, data, associations, jointcalControl)
Definition jointcal.py:768
_iterate_fit(self, associations, fitter, max_steps, name, whatToFit, dataName="", sigmaRelativeTolerance=0, doRankUpdate=True, doLineSearch=False)
Definition jointcal.py:1322
_put_metrics(self, butlerQC, job, outputRefs)
Definition jointcal.py:549
_do_load_refcat_and_fit(self, associations, defaultFilter, center, radius, tract="", match_cut=3.0, reject_bad_fluxes=False, *name="", refObjLoader=None, referenceSelector=None, fit_function=None, epoch=None)
Definition jointcal.py:908
_prep_sky(self, associations, filters)
Definition jointcal.py:818
_compute_proper_motion_epoch(self, ccdImageList)
Definition jointcal.py:887
_make_one_input_data(self, visitRecord, catalog, detectorDict)
Definition jointcal.py:754
_check_star_lists(self, associations, name)
Definition jointcal.py:1047
_fit_photometry(self, associations, dataName=None)
Definition jointcal.py:1099
_put_output(self, butlerQC, outputs, outputRefs, camera, setter)
Definition jointcal.py:565
_load_reference_catalog(self, refObjLoader, referenceSelector, center, radius, filterLabel, applyColorterms=False, epoch=None)
Definition jointcal.py:991
_logChi2AndValidate(self, associations, fit, model, chi2Label, writeChi2Name=None)
Definition jointcal.py:1056
_get_refcat_coordinate_error_override(self, refCat, name)
Definition jointcal.py:840
_load_data(self, inputSourceTableVisit, inputVisitSummary, associations, jointcalControl, camera)
Definition jointcal.py:678
_make_output(self, ccdImageList, model, func)
Definition jointcal.py:1432
runQuantum(self, butlerQC, inputRefs, outputRefs)
Definition jointcal.py:518
_fit_astrometry(self, associations, dataName=None)
Definition jointcal.py:1205
run(self, inputSourceTableVisit, inputVisitSummary, inputCamera, tract=None)
Definition jointcal.py:624
extract_detector_catalog_from_visit_catalog(table, visitCatalog, detectorId, ixxColumns, sourceFluxType, log)
Definition jointcal.py:1535
add_measurement(job, name, value)
Definition jointcal.py:55
writeModel(model, filename, log)
Definition jointcal.py:467
get_sourceTable_visit_columns(inColumns, config, sourceSelector)
Definition jointcal.py:1486
T remove(T... args)