LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
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 
22 import dataclasses
23 import collections
24 import os
25 
26 import astropy.time
27 import numpy as np
28 import astropy.units as u
29 
30 import lsst.geom
31 import lsst.utils
32 import lsst.pex.config as pexConfig
33 import lsst.pipe.base as pipeBase
34 from lsst.afw.image import fluxErrFromABMagErr
35 import lsst.pex.exceptions as pexExceptions
37 import lsst.afw.table
38 import lsst.log
39 from lsst.obs.base import Instrument
40 from lsst.pipe.tasks.colorterms import ColortermLibrary
41 from lsst.verify import Job, Measurement
42 
43 from lsst.meas.algorithms import (LoadIndexedReferenceObjectsTask, ReferenceSourceSelectorTask,
44  ReferenceObjectLoader)
45 from lsst.meas.algorithms.sourceSelector import sourceSelectorRegistry
46 from lsst.utils.timer import timeMethod
47 
48 from .dataIds import PerTractCcdDataIdContainer
49 
50 import lsst.jointcal
51 from lsst.jointcal import MinimizeResult
52 
53 __all__ = ["JointcalConfig", "JointcalRunner", "JointcalTask"]
54 
55 Photometry = collections.namedtuple('Photometry', ('fit', 'model'))
56 Astrometry = collections.namedtuple('Astrometry', ('fit', 'model', 'sky_to_tan_projection'))
57 
58 
59 # TODO: move this to MeasurementSet in lsst.verify per DM-12655.
60 def add_measurement(job, name, value):
61  meas = Measurement(job.metrics[name], value)
62  job.measurements.insert(meas)
63 
64 
65 class JointcalRunner(pipeBase.ButlerInitializedTaskRunner):
66  """Subclass of TaskRunner for jointcalTask (gen2)
67 
68  jointcalTask.runDataRef() takes a number of arguments, one of which is a list of dataRefs
69  extracted from the command line (whereas most CmdLineTasks' runDataRef methods take
70  single dataRef, are are called repeatedly). This class transforms the processed
71  arguments generated by the ArgumentParser into the arguments expected by
72  Jointcal.runDataRef().
73 
74  See pipeBase.TaskRunner for more information.
75  """
76 
77  @staticmethod
78  def getTargetList(parsedCmd, **kwargs):
79  """
80  Return a list of tuples per tract, each containing (dataRefs, kwargs).
81 
82  Jointcal operates on lists of dataRefs simultaneously.
83  """
84  kwargs['butler'] = parsedCmd.butler
85 
86  # organize data IDs by tract
87  refListDict = {}
88  for ref in parsedCmd.id.refList:
89  refListDict.setdefault(ref.dataId["tract"], []).append(ref)
90  # we call runDataRef() once with each tract
91  result = [(refListDict[tract], kwargs) for tract in sorted(refListDict.keys())]
92  return result
93 
94  def __call__(self, args):
95  """
96  Parameters
97  ----------
98  args
99  Arguments for Task.runDataRef()
100 
101  Returns
102  -------
103  pipe.base.Struct
104  if self.doReturnResults is False:
105 
106  - ``exitStatus``: 0 if the task completed successfully, 1 otherwise.
107 
108  if self.doReturnResults is True:
109 
110  - ``result``: the result of calling jointcal.runDataRef()
111  - ``exitStatus``: 0 if the task completed successfully, 1 otherwise.
112  """
113  exitStatus = 0 # exit status for shell
114 
115  # NOTE: cannot call self.makeTask because that assumes args[0] is a single dataRef.
116  dataRefList, kwargs = args
117  butler = kwargs.pop('butler')
118  task = self.TaskClass(config=self.config, log=self.log, butler=butler)
119  result = None
120  try:
121  result = task.runDataRef(dataRefList, **kwargs)
122  exitStatus = result.exitStatus
123  job_path = butler.get('verify_job_filename')
124  result.job.write(job_path[0])
125  except Exception as e: # catch everything, sort it out later.
126  if self.doRaise:
127  raise e
128  else:
129  exitStatus = 1
130  eName = type(e).__name__
131  tract = dataRefList[0].dataId['tract']
132  task.log.fatal("Failed processing tract %s, %s: %s", tract, eName, e)
133 
134  # Put the butler back into kwargs for the other Tasks.
135  kwargs['butler'] = butler
136  if self.doReturnResults:
137  return pipeBase.Struct(result=result, exitStatus=exitStatus)
138  else:
139  return pipeBase.Struct(exitStatus=exitStatus)
140 
141 
142 def lookupStaticCalibrations(datasetType, registry, quantumDataId, collections):
143  """Lookup function that asserts/hopes that a static calibration dataset
144  exists in a particular collection, since this task can't provide a single
145  date/time to use to search for one properly.
146 
147  This is mostly useful for the ``camera`` dataset, in cases where the task's
148  quantum dimensions do *not* include something temporal, like ``exposure``
149  or ``visit``.
150 
151  Parameters
152  ----------
153  datasetType : `lsst.daf.butler.DatasetType`
154  Type of dataset being searched for.
155  registry : `lsst.daf.butler.Registry`
156  Data repository registry to search.
157  quantumDataId : `lsst.daf.butler.DataCoordinate`
158  Data ID of the quantum this camera should match.
159  collections : `Iterable` [ `str` ]
160  Collections that should be searched - but this lookup function works
161  by ignoring this in favor of a more-or-less hard-coded value.
162 
163  Returns
164  -------
165  refs : `Iterator` [ `lsst.daf.butler.DatasetRef` ]
166  Iterator over dataset references; should have only one element.
167 
168  Notes
169  -----
170  This implementation duplicates one in fgcmcal, and is at least quite
171  similar to another in cp_pipe. This duplicate has the most documentation.
172  Fixing this is DM-29661.
173  """
174  instrument = Instrument.fromName(quantumDataId["instrument"], registry)
175  unboundedCollection = instrument.makeUnboundedCalibrationRunName()
176  return registry.queryDatasets(datasetType,
177  dataId=quantumDataId,
178  collections=[unboundedCollection],
179  findFirst=True)
180 
181 
182 def lookupVisitRefCats(datasetType, registry, quantumDataId, collections):
183  """Lookup function that finds all refcats for all visits that overlap a
184  tract, rather than just the refcats that directly overlap the tract.
185 
186  Parameters
187  ----------
188  datasetType : `lsst.daf.butler.DatasetType`
189  Type of dataset being searched for.
190  registry : `lsst.daf.butler.Registry`
191  Data repository registry to search.
192  quantumDataId : `lsst.daf.butler.DataCoordinate`
193  Data ID of the quantum; expected to be something we can use as a
194  constraint to query for overlapping visits.
195  collections : `Iterable` [ `str` ]
196  Collections to search.
197 
198  Returns
199  -------
200  refs : `Iterator` [ `lsst.daf.butler.DatasetRef` ]
201  Iterator over refcat references.
202  """
203  refs = set()
204  # Use .expanded() on the query methods below because we need data IDs with
205  # regions, both in the outer loop over visits (queryDatasets will expand
206  # any data ID we give it, but doing it up-front in bulk is much more
207  # efficient) and in the data IDs of the DatasetRefs this function yields
208  # (because the RefCatLoader relies on them to do some of its own
209  # filtering).
210  for visit_data_id in set(registry.queryDataIds("visit", dataId=quantumDataId).expanded()):
211  refs.update(
212  registry.queryDatasets(
213  datasetType,
214  collections=collections,
215  dataId=visit_data_id,
216  findFirst=True,
217  ).expanded()
218  )
219  yield from refs
220 
221 
222 class JointcalTaskConnections(pipeBase.PipelineTaskConnections,
223  dimensions=("skymap", "tract", "instrument", "physical_filter")):
224  """Middleware input/output connections for jointcal data."""
225  inputCamera = pipeBase.connectionTypes.PrerequisiteInput(
226  doc="The camera instrument that took these observations.",
227  name="camera",
228  storageClass="Camera",
229  dimensions=("instrument",),
230  isCalibration=True,
231  lookupFunction=lookupStaticCalibrations,
232  )
233  inputSourceTableVisit = pipeBase.connectionTypes.Input(
234  doc="Source table in parquet format, per visit",
235  name="sourceTable_visit",
236  storageClass="DataFrame",
237  dimensions=("instrument", "visit"),
238  deferLoad=True,
239  multiple=True,
240  )
241  inputVisitSummary = pipeBase.connectionTypes.Input(
242  doc=("Per-visit consolidated exposure metadata built from calexps. "
243  "These catalogs use detector id for the id and must be sorted for "
244  "fast lookups of a detector."),
245  name="visitSummary",
246  storageClass="ExposureCatalog",
247  dimensions=("instrument", "visit"),
248  deferLoad=True,
249  multiple=True,
250  )
251  astrometryRefCat = pipeBase.connectionTypes.PrerequisiteInput(
252  doc="The astrometry reference catalog to match to loaded input catalog sources.",
253  name="gaia_dr2_20200414",
254  storageClass="SimpleCatalog",
255  dimensions=("skypix",),
256  deferLoad=True,
257  multiple=True,
258  lookupFunction=lookupVisitRefCats,
259  )
260  photometryRefCat = pipeBase.connectionTypes.PrerequisiteInput(
261  doc="The photometry reference catalog to match to loaded input catalog sources.",
262  name="ps1_pv3_3pi_20170110",
263  storageClass="SimpleCatalog",
264  dimensions=("skypix",),
265  deferLoad=True,
266  multiple=True,
267  lookupFunction=lookupVisitRefCats,
268  )
269 
270  outputWcs = pipeBase.connectionTypes.Output(
271  doc=("Per-tract, per-visit world coordinate systems derived from the fitted model."
272  " These catalogs only contain entries for detectors with an output, and use"
273  " the detector id for the catalog id, sorted on id for fast lookups of a detector."),
274  name="jointcalSkyWcsCatalog",
275  storageClass="ExposureCatalog",
276  dimensions=("instrument", "visit", "skymap", "tract"),
277  multiple=True
278  )
279  outputPhotoCalib = pipeBase.connectionTypes.Output(
280  doc=("Per-tract, per-visit photometric calibrations derived from the fitted model."
281  " These catalogs only contain entries for detectors with an output, and use"
282  " the detector id for the catalog id, sorted on id for fast lookups of a detector."),
283  name="jointcalPhotoCalibCatalog",
284  storageClass="ExposureCatalog",
285  dimensions=("instrument", "visit", "skymap", "tract"),
286  multiple=True
287  )
288 
289  # measurements of metrics
290  # The vars() trick used here allows us to set class attributes
291  # programatically. Taken from:
292  # https://stackoverflow.com/questions/2519807/setting-a-class-attribute-with-a-given-name-in-python-while-defining-the-class
293  for name in ("astrometry", "photometry"):
294  vars()[f"{name}_matched_fittedStars"] = pipeBase.connectionTypes.Output(
295  doc=f"The number of cross-matched fittedStars for {name}",
296  name=f"metricvalue_jointcal_{name}_matched_fittedStars",
297  storageClass="MetricValue",
298  dimensions=("skymap", "tract", "instrument", "physical_filter"),
299  )
300  vars()[f"{name}_collected_refStars"] = pipeBase.connectionTypes.Output(
301  doc=f"The number of {name} reference stars drawn from the reference catalog, before matching.",
302  name=f"metricvalue_jointcal_{name}_collected_refStars",
303  storageClass="MetricValue",
304  dimensions=("skymap", "tract", "instrument", "physical_filter"),
305  )
306  vars()[f"{name}_prepared_refStars"] = pipeBase.connectionTypes.Output(
307  doc=f"The number of {name} reference stars matched to fittedStars.",
308  name=f"metricvalue_jointcal_{name}_prepared_refStars",
309  storageClass="MetricValue",
310  dimensions=("skymap", "tract", "instrument", "physical_filter"),
311  )
312  vars()[f"{name}_prepared_fittedStars"] = pipeBase.connectionTypes.Output(
313  doc=f"The number of cross-matched fittedStars after cleanup, for {name}.",
314  name=f"metricvalue_jointcal_{name}_prepared_fittedStars",
315  storageClass="MetricValue",
316  dimensions=("skymap", "tract", "instrument", "physical_filter"),
317  )
318  vars()[f"{name}_prepared_ccdImages"] = pipeBase.connectionTypes.Output(
319  doc=f"The number of ccdImages that will be fit for {name}, after cleanup.",
320  name=f"metricvalue_jointcal_{name}_prepared_ccdImages",
321  storageClass="MetricValue",
322  dimensions=("skymap", "tract", "instrument", "physical_filter"),
323  )
324  vars()[f"{name}_final_chi2"] = pipeBase.connectionTypes.Output(
325  doc=f"The final chi2 of the {name} fit.",
326  name=f"metricvalue_jointcal_{name}_final_chi2",
327  storageClass="MetricValue",
328  dimensions=("skymap", "tract", "instrument", "physical_filter"),
329  )
330  vars()[f"{name}_final_ndof"] = pipeBase.connectionTypes.Output(
331  doc=f"The number of degrees of freedom of the fitted {name}.",
332  name=f"metricvalue_jointcal_{name}_final_ndof",
333  storageClass="MetricValue",
334  dimensions=("skymap", "tract", "instrument", "physical_filter"),
335  )
336 
337  def __init__(self, *, config=None):
338  super().__init__(config=config)
339  # When we are only doing one of astrometry or photometry, we don't
340  # need the reference catalog or produce the outputs for the other.
341  # This informs the middleware of that when the QuantumGraph is
342  # generated, so we don't block on getting something we won't need or
343  # create an expectation that downstream tasks will be able to consume
344  # something we won't produce.
345  if not config.doAstrometry:
346  self.prerequisiteInputs.remove("astrometryRefCat")
347  self.outputs.remove("outputWcs")
348  for key in list(self.outputs):
349  if "metricvalue_jointcal_astrometry" in key:
350  self.outputs.remove(key)
351  if not config.doPhotometry:
352  self.prerequisiteInputs.remove("photometryRefCat")
353  self.outputs.remove("outputPhotoCalib")
354  for key in list(self.outputs):
355  if "metricvalue_jointcal_photometry" in key:
356  self.outputs.remove(key)
357 
358 
359 class JointcalConfig(pipeBase.PipelineTaskConfig,
360  pipelineConnections=JointcalTaskConnections):
361  """Configuration for JointcalTask"""
362 
363  doAstrometry = pexConfig.Field(
364  doc="Fit astrometry and write the fitted result.",
365  dtype=bool,
366  default=True
367  )
368  doPhotometry = pexConfig.Field(
369  doc="Fit photometry and write the fitted result.",
370  dtype=bool,
371  default=True
372  )
373  coaddName = pexConfig.Field(
374  doc="Type of coadd, typically deep or goodSeeing",
375  dtype=str,
376  default="deep"
377  )
378  # TODO DM-29008: Change this to "ApFlux_12_0" before gen2 removal.
379  sourceFluxType = pexConfig.Field(
380  dtype=str,
381  doc="Source flux field to use in source selection and to get fluxes from the catalog.",
382  default='Calib'
383  )
384  positionErrorPedestal = pexConfig.Field(
385  doc="Systematic term to apply to the measured position error (pixels)",
386  dtype=float,
387  default=0.02,
388  )
389  photometryErrorPedestal = pexConfig.Field(
390  doc="Systematic term to apply to the measured error on flux or magnitude as a "
391  "fraction of source flux or magnitude delta (e.g. 0.05 is 5% of flux or +50 millimag).",
392  dtype=float,
393  default=0.0,
394  )
395  # TODO: DM-6885 matchCut should be an geom.Angle
396  matchCut = pexConfig.Field(
397  doc="Matching radius between fitted and reference stars (arcseconds)",
398  dtype=float,
399  default=3.0,
400  )
401  minMeasurements = pexConfig.Field(
402  doc="Minimum number of associated measured stars for a fitted star to be included in the fit",
403  dtype=int,
404  default=2,
405  )
406  minMeasuredStarsPerCcd = pexConfig.Field(
407  doc="Minimum number of measuredStars per ccdImage before printing warnings",
408  dtype=int,
409  default=100,
410  )
411  minRefStarsPerCcd = pexConfig.Field(
412  doc="Minimum number of measuredStars per ccdImage before printing warnings",
413  dtype=int,
414  default=30,
415  )
416  allowLineSearch = pexConfig.Field(
417  doc="Allow a line search during minimization, if it is reasonable for the model"
418  " (models with a significant non-linear component, e.g. constrainedPhotometry).",
419  dtype=bool,
420  default=False
421  )
422  astrometrySimpleOrder = pexConfig.Field(
423  doc="Polynomial order for fitting the simple astrometry model.",
424  dtype=int,
425  default=3,
426  )
427  astrometryChipOrder = pexConfig.Field(
428  doc="Order of the per-chip transform for the constrained astrometry model.",
429  dtype=int,
430  default=1,
431  )
432  astrometryVisitOrder = pexConfig.Field(
433  doc="Order of the per-visit transform for the constrained astrometry model.",
434  dtype=int,
435  default=5,
436  )
437  useInputWcs = pexConfig.Field(
438  doc="Use the input calexp WCSs to initialize a SimpleAstrometryModel.",
439  dtype=bool,
440  default=True,
441  )
442  astrometryModel = pexConfig.ChoiceField(
443  doc="Type of model to fit to astrometry",
444  dtype=str,
445  default="constrained",
446  allowed={"simple": "One polynomial per ccd",
447  "constrained": "One polynomial per ccd, and one polynomial per visit"}
448  )
449  photometryModel = pexConfig.ChoiceField(
450  doc="Type of model to fit to photometry",
451  dtype=str,
452  default="constrainedMagnitude",
453  allowed={"simpleFlux": "One constant zeropoint per ccd and visit, fitting in flux space.",
454  "constrainedFlux": "Constrained zeropoint per ccd, and one polynomial per visit,"
455  " fitting in flux space.",
456  "simpleMagnitude": "One constant zeropoint per ccd and visit,"
457  " fitting in magnitude space.",
458  "constrainedMagnitude": "Constrained zeropoint per ccd, and one polynomial per visit,"
459  " fitting in magnitude space.",
460  }
461  )
462  applyColorTerms = pexConfig.Field(
463  doc="Apply photometric color terms to reference stars?"
464  "Requires that colorterms be set to a ColortermLibrary",
465  dtype=bool,
466  default=False
467  )
468  colorterms = pexConfig.ConfigField(
469  doc="Library of photometric reference catalog name to color term dict.",
470  dtype=ColortermLibrary,
471  )
472  photometryVisitOrder = pexConfig.Field(
473  doc="Order of the per-visit polynomial transform for the constrained photometry model.",
474  dtype=int,
475  default=7,
476  )
477  photometryDoRankUpdate = pexConfig.Field(
478  doc=("Do the rank update step during minimization. "
479  "Skipping this can help deal with models that are too non-linear."),
480  dtype=bool,
481  default=True,
482  )
483  astrometryDoRankUpdate = pexConfig.Field(
484  doc=("Do the rank update step during minimization (should not change the astrometry fit). "
485  "Skipping this can help deal with models that are too non-linear."),
486  dtype=bool,
487  default=True,
488  )
489  outlierRejectSigma = pexConfig.Field(
490  doc="How many sigma to reject outliers at during minimization.",
491  dtype=float,
492  default=5.0,
493  )
494  astrometryOutlierRelativeTolerance = pexConfig.Field(
495  doc=("Convergence tolerance for outlier rejection threshold when fitting astrometry. Iterations will "
496  "stop when the fractional change in the chi2 cut level is below this value. If tolerance is set "
497  "to zero, iterations will continue until there are no more outliers. We suggest a value of 0.002"
498  "as a balance between a shorter minimization runtime and achieving a final fitted model that is"
499  "close to the solution found when removing all outliers."),
500  dtype=float,
501  default=0,
502  )
503  maxPhotometrySteps = pexConfig.Field(
504  doc="Maximum number of minimize iterations to take when fitting photometry.",
505  dtype=int,
506  default=20,
507  )
508  maxAstrometrySteps = pexConfig.Field(
509  doc="Maximum number of minimize iterations to take when fitting astrometry.",
510  dtype=int,
511  default=20,
512  )
513  astrometryRefObjLoader = pexConfig.ConfigurableField(
514  target=LoadIndexedReferenceObjectsTask,
515  doc="Reference object loader for astrometric fit",
516  )
517  photometryRefObjLoader = pexConfig.ConfigurableField(
518  target=LoadIndexedReferenceObjectsTask,
519  doc="Reference object loader for photometric fit",
520  )
521  sourceSelector = sourceSelectorRegistry.makeField(
522  doc="How to select sources for cross-matching",
523  default="astrometry"
524  )
525  astrometryReferenceSelector = pexConfig.ConfigurableField(
526  target=ReferenceSourceSelectorTask,
527  doc="How to down-select the loaded astrometry reference catalog.",
528  )
529  photometryReferenceSelector = pexConfig.ConfigurableField(
530  target=ReferenceSourceSelectorTask,
531  doc="How to down-select the loaded photometry reference catalog.",
532  )
533  astrometryReferenceErr = pexConfig.Field(
534  doc=("Uncertainty on reference catalog coordinates [mas] to use in place of the `coord_*Err` fields. "
535  "If None, then raise an exception if the reference catalog is missing coordinate errors. "
536  "If specified, overrides any existing `coord_*Err` values."),
537  dtype=float,
538  default=None,
539  optional=True
540  )
541 
542  # configs for outputting debug information
543  writeInitMatrix = pexConfig.Field(
544  dtype=bool,
545  doc=("Write the pre/post-initialization Hessian and gradient to text files, for debugging. "
546  "Output files will be written to `config.debugOutputPath` and will "
547  "be of the form 'astrometry_[pre|post]init-TRACT-FILTER-mat.txt'. "
548  "Note that these files are the dense versions of the matrix, and so may be very large."),
549  default=False
550  )
551  writeChi2FilesInitialFinal = pexConfig.Field(
552  dtype=bool,
553  doc=("Write .csv files containing the contributions to chi2 for the initialization and final fit. "
554  "Output files will be written to `config.debugOutputPath` and will "
555  "be of the form `astrometry_[initial|final]_chi2-TRACT-FILTER."),
556  default=False
557  )
558  writeChi2FilesOuterLoop = pexConfig.Field(
559  dtype=bool,
560  doc=("Write .csv files containing the contributions to chi2 for the outer fit loop. "
561  "Output files will be written to `config.debugOutputPath` and will "
562  "be of the form `astrometry_init-NN_chi2-TRACT-FILTER`."),
563  default=False
564  )
565  writeInitialModel = pexConfig.Field(
566  dtype=bool,
567  doc=("Write the pre-initialization model to text files, for debugging. "
568  "Output files will be written to `config.debugOutputPath` and will be "
569  "of the form `initial_astrometry_model-TRACT_FILTER.txt`."),
570  default=False
571  )
572  debugOutputPath = pexConfig.Field(
573  dtype=str,
574  default=".",
575  doc=("Path to write debug output files to. Used by "
576  "`writeInitialModel`, `writeChi2Files*`, `writeInitMatrix`.")
577  )
578  detailedProfile = pexConfig.Field(
579  dtype=bool,
580  default=False,
581  doc="Output separate profiling information for different parts of jointcal, e.g. data read, fitting"
582  )
583 
584  def validate(self):
585  super().validate()
586  if self.doPhotometrydoPhotometry and self.applyColorTermsapplyColorTerms and len(self.colortermscolorterms.data) == 0:
587  msg = "applyColorTerms=True requires the `colorterms` field be set to a ColortermLibrary."
588  raise pexConfig.FieldValidationError(JointcalConfig.colorterms, self, msg)
589  if self.doAstrometrydoAstrometry and not self.doPhotometrydoPhotometry and self.applyColorTermsapplyColorTerms:
590  msg = ("Only doing astrometry, but Colorterms are not applied for astrometry;"
591  "applyColorTerms=True will be ignored.")
592  lsst.log.warning(msg)
593 
594  def setDefaults(self):
595  # Use science source selector which can filter on extendedness, SNR, and whether blended
596  self.sourceSelectorsourceSelector.name = 'science'
597  # Use only stars because aperture fluxes of galaxies are biased and depend on seeing
598  self.sourceSelectorsourceSelector['science'].doUnresolved = True
599  # with dependable signal to noise ratio.
600  self.sourceSelectorsourceSelector['science'].doSignalToNoise = True
601  # Min SNR must be > 0 because jointcal cannot handle negative fluxes,
602  # and S/N > 10 to use sources that are not too faint, and thus better measured.
603  self.sourceSelectorsourceSelector['science'].signalToNoise.minimum = 10.
604  # Base SNR on CalibFlux because that is the flux jointcal that fits and must be positive
605  fluxField = f"slot_{self.sourceFluxType}Flux_instFlux"
606  self.sourceSelectorsourceSelector['science'].signalToNoise.fluxField = fluxField
607  self.sourceSelectorsourceSelector['science'].signalToNoise.errField = fluxField + "Err"
608  # Do not trust blended sources' aperture fluxes which also depend on seeing
609  self.sourceSelectorsourceSelector['science'].doIsolated = True
610  # Do not trust either flux or centroid measurements with flags,
611  # chosen from the usual QA flags for stars)
612  self.sourceSelectorsourceSelector['science'].doFlags = True
613  badFlags = ['base_PixelFlags_flag_edge', 'base_PixelFlags_flag_saturated',
614  'base_PixelFlags_flag_interpolatedCenter', 'base_SdssCentroid_flag',
615  'base_PsfFlux_flag', 'base_PixelFlags_flag_suspectCenter']
616  self.sourceSelectorsourceSelector['science'].flags.bad = badFlags
617 
618  # Default to Gaia-DR2 (with proper motions) for astrometry and
619  # PS1-DR1 for photometry, with a reasonable initial filterMap.
620  self.astrometryRefObjLoaderastrometryRefObjLoader.ref_dataset_name = "gaia_dr2_20200414"
621  self.astrometryRefObjLoaderastrometryRefObjLoader.requireProperMotion = True
622  self.astrometryRefObjLoaderastrometryRefObjLoader.anyFilterMapsToThis = 'phot_g_mean'
623  self.photometryRefObjLoaderphotometryRefObjLoader.ref_dataset_name = "ps1_pv3_3pi_20170110"
624 
625 
626 def writeModel(model, filename, log):
627  """Write model to outfile."""
628  with open(filename, "w") as file:
629  file.write(repr(model))
630  log.info("Wrote %s to file: %s", model, filename)
631 
632 
633 @dataclasses.dataclass
635  """The input data jointcal needs for each detector/visit."""
636  visit: int
637  """The visit identifier of this exposure."""
639  """The catalog derived from this exposure."""
640  visitInfo: lsst.afw.image.VisitInfo
641  """The VisitInfo of this exposure."""
643  """The detector of this exposure."""
644  photoCalib: lsst.afw.image.PhotoCalib
645  """The photometric calibration of this exposure."""
647  """The WCS of this exposure."""
648  bbox: lsst.geom.Box2I
649  """The bounding box of this exposure."""
651  """The filter of this exposure."""
652 
653 
654 class JointcalTask(pipeBase.PipelineTask, pipeBase.CmdLineTask):
655  """Astrometricly and photometricly calibrate across multiple visits of the
656  same field.
657 
658  Parameters
659  ----------
660  butler : `lsst.daf.persistence.Butler`
661  The butler is passed to the refObjLoader constructor in case it is
662  needed. Ignored if the refObjLoader argument provides a loader directly.
663  Used to initialize the astrometry and photometry refObjLoaders.
664  initInputs : `dict`, optional
665  Dictionary used to initialize PipelineTasks (empty for jointcal).
666  """
667 
668  ConfigClass = JointcalConfig
669  RunnerClass = JointcalRunner
670  _DefaultName = "jointcal"
671 
672  def __init__(self, butler=None, initInputs=None, **kwargs):
673  super().__init__(**kwargs)
674  self.makeSubtask("sourceSelector")
675  if self.config.doAstrometry:
676  if initInputs is None:
677  # gen3 middleware does refcat things internally (and will not have a butler here)
678  self.makeSubtask('astrometryRefObjLoader', butler=butler)
679  self.makeSubtask("astrometryReferenceSelector")
680  else:
681  self.astrometryRefObjLoaderastrometryRefObjLoader = None
682  if self.config.doPhotometry:
683  if initInputs is None:
684  # gen3 middleware does refcat things internally (and will not have a butler here)
685  self.makeSubtask('photometryRefObjLoader', butler=butler)
686  self.makeSubtask("photometryReferenceSelector")
687  else:
688  self.photometryRefObjLoaderphotometryRefObjLoader = None
689 
690  # To hold various computed metrics for use by tests
691  self.jobjob = Job.load_metrics_package(subset='jointcal')
692 
693  def runQuantum(self, butlerQC, inputRefs, outputRefs):
694  # We override runQuantum to set up the refObjLoaders and write the
695  # outputs to the correct refs.
696  inputs = butlerQC.get(inputRefs)
697  # We want the tract number for writing debug files
698  tract = butlerQC.quantum.dataId['tract']
699  if self.config.doAstrometry:
700  self.astrometryRefObjLoaderastrometryRefObjLoader = ReferenceObjectLoader(
701  dataIds=[ref.datasetRef.dataId
702  for ref in inputRefs.astrometryRefCat],
703  refCats=inputs.pop('astrometryRefCat'),
704  config=self.config.astrometryRefObjLoader,
705  log=self.log)
706  if self.config.doPhotometry:
707  self.photometryRefObjLoaderphotometryRefObjLoader = ReferenceObjectLoader(
708  dataIds=[ref.datasetRef.dataId
709  for ref in inputRefs.photometryRefCat],
710  refCats=inputs.pop('photometryRefCat'),
711  config=self.config.photometryRefObjLoader,
712  log=self.log)
713  outputs = self.runrun(**inputs, tract=tract)
714  self._put_metrics_put_metrics(butlerQC, outputs.job, outputRefs)
715  if self.config.doAstrometry:
716  self._put_output_put_output(butlerQC, outputs.outputWcs, outputRefs.outputWcs,
717  inputs['inputCamera'], "setWcs")
718  if self.config.doPhotometry:
719  self._put_output_put_output(butlerQC, outputs.outputPhotoCalib, outputRefs.outputPhotoCalib,
720  inputs['inputCamera'], "setPhotoCalib")
721 
722  def _put_metrics(self, butlerQC, job, outputRefs):
723  """Persist all measured metrics stored in a job.
724 
725  Parameters
726  ----------
727  butlerQC : `lsst.pipe.base.ButlerQuantumContext`
728  A butler which is specialized to operate in the context of a
729  `lsst.daf.butler.Quantum`; This is the input to `runQuantum`.
730  job : `lsst.verify.job.Job`
731  Measurements of metrics to persist.
732  outputRefs : `list` [`lsst.pipe.base.connectionTypes.OutputQuantizedConnection`]
733  The DatasetRefs to persist the data to.
734  """
735  for key in job.measurements.keys():
736  butlerQC.put(job.measurements[key], getattr(outputRefs, key.fqn.replace('jointcal.', '')))
737 
738  def _put_output(self, butlerQC, outputs, outputRefs, camera, setter):
739  """Persist the output datasets to their appropriate datarefs.
740 
741  Parameters
742  ----------
743  butlerQC : `lsst.pipe.base.ButlerQuantumContext`
744  A butler which is specialized to operate in the context of a
745  `lsst.daf.butler.Quantum`; This is the input to `runQuantum`.
746  outputs : `dict` [`tuple`, `lsst.afw.geom.SkyWcs`] or
747  `dict` [`tuple, `lsst.afw.image.PhotoCalib`]
748  The fitted objects to persist.
749  outputRefs : `list` [`lsst.pipe.base.connectionTypes.OutputQuantizedConnection`]
750  The DatasetRefs to persist the data to.
751  camera : `lsst.afw.cameraGeom.Camera`
752  The camera for this instrument, to get detector ids from.
753  setter : `str`
754  The method to call on the ExposureCatalog to set each output.
755  """
757  schema.addField('visit', type='I', doc='Visit number')
758 
759  def new_catalog(visit, size):
760  """Return an catalog ready to be filled with appropriate output."""
761  catalog = lsst.afw.table.ExposureCatalog(schema)
762  catalog.resize(size)
763  catalog['visit'] = visit
764  metadata = lsst.daf.base.PropertyList()
765  metadata.add("COMMENT", "Catalog id is detector id, sorted.")
766  metadata.add("COMMENT", "Only detectors with data have entries.")
767  return catalog
768 
769  # count how many detectors have output for each visit
770  detectors_per_visit = collections.defaultdict(int)
771  for key in outputs:
772  # key is (visit, detector_id), and we only need visit here
773  detectors_per_visit[key[0]] += 1
774 
775  for ref in outputRefs:
776  visit = ref.dataId['visit']
777  catalog = new_catalog(visit, detectors_per_visit[visit])
778  # Iterate over every detector and skip the ones we don't have output for.
779  i = 0
780  for detector in camera:
781  detectorId = detector.getId()
782  key = (ref.dataId['visit'], detectorId)
783  if key not in outputs:
784  # skip detectors we don't have output for
785  self.log.debug("No %s output for detector %s in visit %s",
786  setter[3:], detectorId, visit)
787  continue
788 
789  catalog[i].setId(detectorId)
790  getattr(catalog[i], setter)(outputs[key])
791  i += 1
792 
793  catalog.sort() # ensure that the detectors are in sorted order, for fast lookups
794  butlerQC.put(catalog, ref)
795  self.log.info("Wrote %s detectors to %s", i, ref)
796 
797  def run(self, inputSourceTableVisit, inputVisitSummary, inputCamera, tract=None):
798  # Docstring inherited.
799 
800  # We take values out of the Parquet table, and put them in "flux_",
801  # and the config.sourceFluxType field is used during that extraction,
802  # so just use "flux" here.
803  sourceFluxField = "flux"
804  jointcalControl = lsst.jointcal.JointcalControl(sourceFluxField)
805  associations = lsst.jointcal.Associations()
806  self.focalPlaneBBoxfocalPlaneBBox = inputCamera.getFpBBox()
807  oldWcsList, bands = self._load_data_load_data(inputSourceTableVisit,
808  inputVisitSummary,
809  associations,
810  jointcalControl,
811  inputCamera)
812 
813  boundingCircle, center, radius, defaultFilter, epoch = self._prep_sky_prep_sky(associations, bands)
814 
815  if self.config.doAstrometry:
816  astrometry = self._do_load_refcat_and_fit_do_load_refcat_and_fit(associations, defaultFilter, center, radius,
817  name="astrometry",
818  refObjLoader=self.astrometryRefObjLoaderastrometryRefObjLoader,
819  referenceSelector=self.astrometryReferenceSelector,
820  fit_function=self._fit_astrometry_fit_astrometry,
821  tract=tract,
822  epoch=epoch)
823  astrometry_output = self._make_output_make_output(associations.getCcdImageList(),
824  astrometry.model,
825  "makeSkyWcs")
826  else:
827  astrometry_output = None
828 
829  if self.config.doPhotometry:
830  photometry = self._do_load_refcat_and_fit_do_load_refcat_and_fit(associations, defaultFilter, center, radius,
831  name="photometry",
832  refObjLoader=self.photometryRefObjLoaderphotometryRefObjLoader,
833  referenceSelector=self.photometryReferenceSelector,
834  fit_function=self._fit_photometry_fit_photometry,
835  tract=tract,
836  epoch=epoch,
837  reject_bad_fluxes=True)
838  photometry_output = self._make_output_make_output(associations.getCcdImageList(),
839  photometry.model,
840  "toPhotoCalib")
841  else:
842  photometry_output = None
843 
844  return pipeBase.Struct(outputWcs=astrometry_output,
845  outputPhotoCalib=photometry_output,
846  job=self.jobjob,
847  astrometryRefObjLoader=self.astrometryRefObjLoaderastrometryRefObjLoader,
848  photometryRefObjLoader=self.photometryRefObjLoaderphotometryRefObjLoader)
849 
850  def _make_schema_table(self):
851  """Return an afw SourceTable to use as a base for creating the
852  SourceCatalog to insert values from the dataFrame into.
853 
854  Returns
855  -------
856  table : `lsst.afw.table.SourceTable`
857  Table with schema and slots to use to make SourceCatalogs.
858  """
860  schema.addField("centroid_x", "D")
861  schema.addField("centroid_y", "D")
862  schema.addField("centroid_xErr", "F")
863  schema.addField("centroid_yErr", "F")
864  schema.addField("shape_xx", "D")
865  schema.addField("shape_yy", "D")
866  schema.addField("shape_xy", "D")
867  schema.addField("flux_instFlux", "D")
868  schema.addField("flux_instFluxErr", "D")
869  table = lsst.afw.table.SourceTable.make(schema)
870  table.defineCentroid("centroid")
871  table.defineShape("shape")
872  return table
873 
874  def _extract_detector_catalog_from_visit_catalog(self, table, visitCatalog, detectorId,
875  detectorColumn, ixxColumns):
876  """Return an afw SourceCatalog extracted from a visit-level dataframe,
877  limited to just one detector.
878 
879  Parameters
880  ----------
881  table : `lsst.afw.table.SourceTable`
882  Table factory to use to make the SourceCatalog that will be
883  populated with data from ``visitCatalog``.
884  visitCatalog : `pandas.DataFrame`
885  DataFrame to extract a detector catalog from.
886  detectorId : `int`
887  Numeric id of the detector to extract from ``visitCatalog``.
888  detectorColumn : `str`
889  Name of the detector column in the catalog.
890  ixxColumns : `list` [`str`]
891  Names of the ixx/iyy/ixy columns in the catalog.
892 
893  Returns
894  -------
895  catalog : `lsst.afw.table.SourceCatalog`
896  Detector-level catalog extracted from ``visitCatalog``.
897  """
898  # map from dataFrame column to afw table column
899  mapping = {'x': 'centroid_x',
900  'y': 'centroid_y',
901  'xErr': 'centroid_xErr',
902  'yErr': 'centroid_yErr',
903  ixxColumns[0]: 'shape_xx',
904  ixxColumns[1]: 'shape_yy',
905  ixxColumns[2]: 'shape_xy',
906  f'{self.config.sourceFluxType}_instFlux': 'flux_instFlux',
907  f'{self.config.sourceFluxType}_instFluxErr': 'flux_instFluxErr',
908  }
909 
910  catalog = lsst.afw.table.SourceCatalog(table)
911  matched = visitCatalog[detectorColumn] == detectorId
912  catalog.resize(sum(matched))
913  view = visitCatalog.loc[matched]
914  catalog['id'] = view.index
915  for dfCol, afwCol in mapping.items():
916  catalog[afwCol] = view[dfCol]
917 
918  self.log.debug("%d sources selected in visit %d detector %d",
919  len(catalog),
920  view['visit'].iloc[0], # all visits in this catalog are the same, so take the first
921  detectorId)
922  return catalog
923 
924  def _load_data(self, inputSourceTableVisit, inputVisitSummary, associations,
925  jointcalControl, camera):
926  """Read the data that jointcal needs to run. (Gen3 version)
927 
928  Modifies ``associations`` in-place with the loaded data.
929 
930  Parameters
931  ----------
932  inputSourceTableVisit : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
933  References to visit-level DataFrames to load the catalog data from.
934  inputVisitSummary : `list` [`lsst.daf.butler.DeferredDatasetHandle`]
935  Visit-level exposure summary catalog with metadata.
936  associations : `lsst.jointcal.Associations`
937  Object to add the loaded data to by constructing new CcdImages.
938  jointcalControl : `jointcal.JointcalControl`
939  Control object for C++ associations management.
940  camera : `lsst.afw.cameraGeom.Camera`
941  Camera object for detector geometry.
942 
943  Returns
944  -------
945  oldWcsList: `list` [`lsst.afw.geom.SkyWcs`]
946  The original WCS of the input data, to aid in writing tests.
947  bands : `list` [`str`]
948  The filter bands of each input dataset.
949  """
950  oldWcsList = []
951  filters = []
952  load_cat_prof_file = 'jointcal_load_data.prof' if self.config.detailedProfile else ''
953  with pipeBase.cmdLineTask.profile(load_cat_prof_file):
954  table = self._make_schema_table_make_schema_table() # every detector catalog has the same layout
955  # No guarantee that the input is in the same order of visits, so we have to map one of them.
956  catalogMap = {ref.dataId['visit']: i for i, ref in enumerate(inputSourceTableVisit)}
957  detectorDict = {detector.getId(): detector for detector in camera}
958 
959  columns = None
960 
961  for visitSummaryRef in inputVisitSummary:
962  visitSummary = visitSummaryRef.get()
963 
964  dataRef = inputSourceTableVisit[catalogMap[visitSummaryRef.dataId['visit']]]
965  if columns is None:
966  inColumns = dataRef.get(component='columns')
967  columns, detColumn, ixxColumns = self._get_sourceTable_visit_columns_get_sourceTable_visit_columns(inColumns)
968  visitCatalog = dataRef.get(parameters={'columns': columns})
969 
970  selected = self.sourceSelector.run(visitCatalog)
971 
972  # Build a CcdImage for each detector in this visit.
973  detectors = {id: index for index, id in enumerate(visitSummary['id'])}
974  for id, index in detectors.items():
975  catalog = self._extract_detector_catalog_from_visit_catalog_extract_detector_catalog_from_visit_catalog(table, selected.sourceCat, id,
976  detColumn, ixxColumns)
977  data = self._make_one_input_data_make_one_input_data(visitSummary[index], catalog, detectorDict)
978  result = self._build_ccdImage_build_ccdImage(data, associations, jointcalControl)
979  if result is not None:
980  oldWcsList.append(result.wcs)
981  # A visit has only one band, so we can just use the first.
982  filters.append(data.filter)
983  if len(filters) == 0:
984  raise RuntimeError("No data to process: did source selector remove all sources?")
985  filters = collections.Counter(filters)
986 
987  return oldWcsList, filters
988 
989  def _make_one_input_data(self, visitRecord, catalog, detectorDict):
990  """Return a data structure for this detector+visit."""
991  return JointcalInputData(visit=visitRecord['visit'],
992  catalog=catalog,
993  visitInfo=visitRecord.getVisitInfo(),
994  detector=detectorDict[visitRecord.getId()],
995  photoCalib=visitRecord.getPhotoCalib(),
996  wcs=visitRecord.getWcs(),
997  bbox=visitRecord.getBBox(),
998  # ExposureRecord doesn't have a FilterLabel yet,
999  # so we have to make one.
1000  filter=lsst.afw.image.FilterLabel(band=visitRecord['band'],
1001  physical=visitRecord['physical_filter']))
1002 
1003  def _get_sourceTable_visit_columns(self, inColumns):
1004  """
1005  Get the sourceTable_visit columns from the config.
1006 
1007  Parameters
1008  ----------
1009  inColumns : `list`
1010  List of columns available in the sourceTable_visit
1011 
1012  Returns
1013  -------
1014  columns : `list`
1015  List of columns to read from sourceTable_visit.
1016  detectorColumn : `str`
1017  Name of the detector column.
1018  ixxColumns : `list`
1019  Name of the ixx/iyy/ixy columns.
1020  """
1021  if 'detector' in inColumns:
1022  # Default name for Gen3.
1023  detectorColumn = 'detector'
1024  else:
1025  # Default name for Gen2 and Gen2 conversions.
1026  detectorColumn = 'ccd'
1027 
1028  columns = ['visit', detectorColumn,
1029  'sourceId', 'x', 'xErr', 'y', 'yErr',
1030  self.config.sourceFluxType + '_instFlux', self.config.sourceFluxType + '_instFluxErr']
1031 
1032  if 'ixx' in inColumns:
1033  # New columns post-DM-31825
1034  ixxColumns = ['ixx', 'iyy', 'ixy']
1035  else:
1036  # Old columns pre-DM-31825
1037  ixxColumns = ['Ixx', 'Iyy', 'Ixy']
1038  columns.extend(ixxColumns)
1039 
1040  if self.sourceSelector.config.doFlags:
1041  columns.extend(self.sourceSelector.config.flags.bad)
1042  if self.sourceSelector.config.doUnresolved:
1043  columns.append(self.sourceSelector.config.unresolved.name)
1044  if self.sourceSelector.config.doIsolated:
1045  columns.append(self.sourceSelector.config.isolated.parentName)
1046  columns.append(self.sourceSelector.config.isolated.nChildName)
1047 
1048  return columns, detectorColumn, ixxColumns
1049 
1050  # We don't currently need to persist the metadata.
1051  # If we do in the future, we will have to add appropriate dataset templates
1052  # to each obs package (the metadata template should look like `jointcal_wcs`).
1053  def _getMetadataName(self):
1054  return None
1055 
1056  @classmethod
1057  def _makeArgumentParser(cls):
1058  """Create an argument parser"""
1059  parser = pipeBase.ArgumentParser(name=cls._DefaultName_DefaultName)
1060  parser.add_id_argument("--id", "calexp", help="data ID, e.g. --id visit=6789 ccd=0..9",
1061  ContainerClass=PerTractCcdDataIdContainer)
1062  return parser
1063 
1064  def _build_ccdImage(self, data, associations, jointcalControl):
1065  """
1066  Extract the necessary things from this catalog+metadata to add a new
1067  ccdImage.
1068 
1069  Parameters
1070  ----------
1071  data : `JointcalInputData`
1072  The loaded input data.
1073  associations : `lsst.jointcal.Associations`
1074  Object to add the info to, to construct a new CcdImage
1075  jointcalControl : `jointcal.JointcalControl`
1076  Control object for associations management
1077 
1078  Returns
1079  ------
1080  namedtuple or `None`
1081  ``wcs``
1082  The TAN WCS of this image, read from the calexp
1083  (`lsst.afw.geom.SkyWcs`).
1084  ``key``
1085  A key to identify this dataRef by its visit and ccd ids
1086  (`namedtuple`).
1087  `None`
1088  if there are no sources in the loaded catalog.
1089  """
1090  if len(data.catalog) == 0:
1091  self.log.warning("No sources selected in visit %s ccd %s", data.visit, data.detector.getId())
1092  return None
1093 
1094  associations.createCcdImage(data.catalog,
1095  data.wcs,
1096  data.visitInfo,
1097  data.bbox,
1098  data.filter.physicalLabel,
1099  data.photoCalib,
1100  data.detector,
1101  data.visit,
1102  data.detector.getId(),
1103  jointcalControl)
1104 
1105  Result = collections.namedtuple('Result_from_build_CcdImage', ('wcs', 'key'))
1106  Key = collections.namedtuple('Key', ('visit', 'ccd'))
1107  return Result(data.wcs, Key(data.visit, data.detector.getId()))
1108 
1109  def _readDataId(self, butler, dataId):
1110  """Read all of the data for one dataId from the butler. (gen2 version)"""
1111  # Not all instruments have `visit` in their dataIds.
1112  if "visit" in dataId.keys():
1113  visit = dataId["visit"]
1114  else:
1115  visit = butler.getButler().queryMetadata("calexp", ("visit"), butler.dataId)[0]
1116  detector = butler.get('calexp_detector', dataId=dataId)
1117 
1118  catalog = butler.get('src',
1119  flags=lsst.afw.table.SOURCE_IO_NO_FOOTPRINTS,
1120  dataId=dataId)
1121  goodSrc = self.sourceSelector.run(catalog)
1122  self.log.debug("%d sources selected in visit %d detector %d",
1123  len(goodSrc.sourceCat),
1124  visit,
1125  detector.getId())
1126  return JointcalInputData(visit=visit,
1127  catalog=goodSrc.sourceCat,
1128  visitInfo=butler.get('calexp_visitInfo', dataId=dataId),
1129  detector=detector,
1130  photoCalib=butler.get('calexp_photoCalib', dataId=dataId),
1131  wcs=butler.get('calexp_wcs', dataId=dataId),
1132  bbox=butler.get('calexp_bbox', dataId=dataId),
1133  filter=butler.get('calexp_filterLabel', dataId=dataId))
1134 
1135  def loadData(self, dataRefs, associations, jointcalControl):
1136  """Read the data that jointcal needs to run. (Gen2 version)"""
1137  visit_ccd_to_dataRef = {}
1138  oldWcsList = []
1139  filters = []
1140  load_cat_prof_file = 'jointcal_loadData.prof' if self.config.detailedProfile else ''
1141  with pipeBase.cmdLineTask.profile(load_cat_prof_file):
1142  # Need the bounding-box of the focal plane (the same for all visits) for photometry visit models
1143  camera = dataRefs[0].get('camera', immediate=True)
1144  self.focalPlaneBBoxfocalPlaneBBox = camera.getFpBBox()
1145  for dataRef in dataRefs:
1146  data = self._readDataId_readDataId(dataRef.getButler(), dataRef.dataId)
1147  result = self._build_ccdImage_build_ccdImage(data, associations, jointcalControl)
1148  if result is None:
1149  continue
1150  oldWcsList.append(result.wcs)
1151  visit_ccd_to_dataRef[result.key] = dataRef
1152  filters.append(data.filter)
1153  if len(filters) == 0:
1154  raise RuntimeError("No data to process: did source selector remove all sources?")
1155  filters = collections.Counter(filters)
1156 
1157  return oldWcsList, filters, visit_ccd_to_dataRef
1158 
1159  def _getDebugPath(self, filename):
1160  """Constructs a path to filename using the configured debug path.
1161  """
1162  return os.path.join(self.config.debugOutputPath, filename)
1163 
1164  def _prep_sky(self, associations, filters):
1165  """Prepare on-sky and other data that must be computed after data has
1166  been read.
1167  """
1168  associations.computeCommonTangentPoint()
1169 
1170  boundingCircle = associations.computeBoundingCircle()
1171  center = lsst.geom.SpherePoint(boundingCircle.getCenter())
1172  radius = lsst.geom.Angle(boundingCircle.getOpeningAngle().asRadians(), lsst.geom.radians)
1173 
1174  self.log.info(f"Data has center={center} with radius={radius.asDegrees()} degrees.")
1175 
1176  # Determine a default filter band associated with the catalog. See DM-9093
1177  defaultFilter = filters.most_common(1)[0][0]
1178  self.log.debug("Using '%s' filter for reference flux", defaultFilter.physicalLabel)
1179 
1180  # compute and set the reference epoch of the observations, for proper motion corrections
1181  epoch = self._compute_proper_motion_epoch_compute_proper_motion_epoch(associations.getCcdImageList())
1182  associations.setEpoch(epoch.jyear)
1183 
1184  return boundingCircle, center, radius, defaultFilter, epoch
1185 
1186  @timeMethod
1187  def runDataRef(self, dataRefs):
1188  """
1189  Jointly calibrate the astrometry and photometry across a set of images.
1190 
1191  NOTE: this is for gen2 middleware only.
1192 
1193  Parameters
1194  ----------
1195  dataRefs : `list` of `lsst.daf.persistence.ButlerDataRef`
1196  List of data references to the exposures to be fit.
1197 
1198  Returns
1199  -------
1200  result : `lsst.pipe.base.Struct`
1201  Struct of metadata from the fit, containing:
1202 
1203  ``dataRefs``
1204  The provided data references that were fit (with updated WCSs)
1205  ``oldWcsList``
1206  The original WCS from each dataRef
1207  ``metrics``
1208  Dictionary of internally-computed metrics for testing/validation.
1209  """
1210  if len(dataRefs) == 0:
1211  raise ValueError('Need a non-empty list of data references!')
1212 
1213  exitStatus = 0 # exit status for shell
1214 
1215  sourceFluxField = "slot_%sFlux" % (self.config.sourceFluxType,)
1216  jointcalControl = lsst.jointcal.JointcalControl(sourceFluxField)
1217  associations = lsst.jointcal.Associations()
1218 
1219  oldWcsList, filters, visit_ccd_to_dataRef = self.loadDataloadData(dataRefs,
1220  associations,
1221  jointcalControl)
1222 
1223  boundingCircle, center, radius, defaultFilter, epoch = self._prep_sky_prep_sky(associations, filters)
1224 
1225  tract = dataRefs[0].dataId['tract']
1226 
1227  if self.config.doAstrometry:
1228  astrometry = self._do_load_refcat_and_fit_do_load_refcat_and_fit(associations, defaultFilter, center, radius,
1229  name="astrometry",
1230  refObjLoader=self.astrometryRefObjLoaderastrometryRefObjLoader,
1231  referenceSelector=self.astrometryReferenceSelector,
1232  fit_function=self._fit_astrometry_fit_astrometry,
1233  tract=tract,
1234  epoch=epoch)
1235  self._write_astrometry_results_write_astrometry_results(associations, astrometry.model, visit_ccd_to_dataRef)
1236  else:
1237  astrometry = Astrometry(None, None, None)
1238 
1239  if self.config.doPhotometry:
1240  photometry = self._do_load_refcat_and_fit_do_load_refcat_and_fit(associations, defaultFilter, center, radius,
1241  name="photometry",
1242  refObjLoader=self.photometryRefObjLoaderphotometryRefObjLoader,
1243  referenceSelector=self.photometryReferenceSelector,
1244  fit_function=self._fit_photometry_fit_photometry,
1245  tract=tract,
1246  epoch=epoch,
1247  reject_bad_fluxes=True)
1248  self._write_photometry_results_write_photometry_results(associations, photometry.model, visit_ccd_to_dataRef)
1249  else:
1250  photometry = Photometry(None, None)
1251 
1252  return pipeBase.Struct(dataRefs=dataRefs,
1253  oldWcsList=oldWcsList,
1254  job=self.jobjob,
1255  astrometryRefObjLoader=self.astrometryRefObjLoaderastrometryRefObjLoader,
1256  photometryRefObjLoader=self.photometryRefObjLoaderphotometryRefObjLoader,
1257  defaultFilter=defaultFilter,
1258  epoch=epoch,
1259  exitStatus=exitStatus)
1260 
1261  def _get_refcat_coordinate_error_override(self, refCat, name):
1262  """Check whether we should override the refcat coordinate errors, and
1263  return the overridden error if necessary.
1264 
1265  Parameters
1266  ----------
1267  refCat : `lsst.afw.table.SimpleCatalog`
1268  The reference catalog to check for a ``coord_raErr`` field.
1269  name : `str`
1270  Whether we are doing "astrometry" or "photometry".
1271 
1272  Returns
1273  -------
1274  refCoordErr : `float`
1275  The refcat coordinate error to use, or NaN if we are not overriding
1276  those fields.
1277 
1278  Raises
1279  ------
1280  lsst.pex.config.FieldValidationError
1281  Raised if the refcat does not contain coordinate errors and
1282  ``config.astrometryReferenceErr`` is not set.
1283  """
1284  # This value doesn't matter for photometry, so just set something to
1285  # keep old refcats from causing problems.
1286  if name.lower() == "photometry":
1287  if 'coord_raErr' not in refCat.schema:
1288  return 100
1289  else:
1290  return float('nan')
1291 
1292  if self.config.astrometryReferenceErr is None and 'coord_raErr' not in refCat.schema:
1293  msg = ("Reference catalog does not contain coordinate errors, "
1294  "and config.astrometryReferenceErr not supplied.")
1295  raise pexConfig.FieldValidationError(JointcalConfig.astrometryReferenceErr,
1296  self.config,
1297  msg)
1298 
1299  if self.config.astrometryReferenceErr is not None and 'coord_raErr' in refCat.schema:
1300  self.log.warning("Overriding reference catalog coordinate errors with %f/coordinate [mas]",
1301  self.config.astrometryReferenceErr)
1302 
1303  if self.config.astrometryReferenceErr is None:
1304  return float('nan')
1305  else:
1306  return self.config.astrometryReferenceErr
1307 
1308  def _compute_proper_motion_epoch(self, ccdImageList):
1309  """Return the proper motion correction epoch of the provided images.
1310 
1311  Parameters
1312  ----------
1313  ccdImageList : `list` [`lsst.jointcal.CcdImage`]
1314  The images to compute the appropriate epoch for.
1315 
1316  Returns
1317  -------
1318  epoch : `astropy.time.Time`
1319  The date to use for proper motion corrections.
1320  """
1321  return astropy.time.Time(np.mean([ccdImage.getEpoch() for ccdImage in ccdImageList]),
1322  format="jyear",
1323  scale="tai")
1324 
1325  def _do_load_refcat_and_fit(self, associations, defaultFilter, center, radius,
1326  tract="", match_cut=3.0,
1327  reject_bad_fluxes=False, *,
1328  name="", refObjLoader=None, referenceSelector=None,
1329  fit_function=None, epoch=None):
1330  """Load reference catalog, perform the fit, and return the result.
1331 
1332  Parameters
1333  ----------
1334  associations : `lsst.jointcal.Associations`
1335  The star/reference star associations to fit.
1336  defaultFilter : `lsst.afw.image.FilterLabel`
1337  filter to load from reference catalog.
1338  center : `lsst.geom.SpherePoint`
1339  ICRS center of field to load from reference catalog.
1340  radius : `lsst.geom.Angle`
1341  On-sky radius to load from reference catalog.
1342  name : `str`
1343  Name of thing being fit: "astrometry" or "photometry".
1344  refObjLoader : `lsst.meas.algorithms.LoadReferenceObjectsTask`
1345  Reference object loader to use to load a reference catalog.
1346  referenceSelector : `lsst.meas.algorithms.ReferenceSourceSelectorTask`
1347  Selector to use to pick objects from the loaded reference catalog.
1348  fit_function : callable
1349  Function to call to perform fit (takes Associations object).
1350  tract : `str`, optional
1351  Name of tract currently being fit.
1352  match_cut : `float`, optional
1353  Radius in arcseconds to find cross-catalog matches to during
1354  associations.associateCatalogs.
1355  reject_bad_fluxes : `bool`, optional
1356  Reject refCat sources with NaN/inf flux or NaN/0 fluxErr.
1357  epoch : `astropy.time.Time`, optional
1358  Epoch to which to correct refcat proper motion and parallax,
1359  or `None` to not apply such corrections.
1360 
1361  Returns
1362  -------
1363  result : `Photometry` or `Astrometry`
1364  Result of `fit_function()`
1365  """
1366  self.log.info("====== Now processing %s...", name)
1367  # TODO: this should not print "trying to invert a singular transformation:"
1368  # if it does that, something's not right about the WCS...
1369  associations.associateCatalogs(match_cut)
1370  add_measurement(self.jobjob, 'jointcal.%s_matched_fittedStars' % name,
1371  associations.fittedStarListSize())
1372 
1373  applyColorterms = False if name.lower() == "astrometry" else self.config.applyColorTerms
1374  refCat, fluxField = self._load_reference_catalog_load_reference_catalog(refObjLoader, referenceSelector,
1375  center, radius, defaultFilter,
1376  applyColorterms=applyColorterms,
1377  epoch=epoch)
1378  refCoordErr = self._get_refcat_coordinate_error_override_get_refcat_coordinate_error_override(refCat, name)
1379 
1380  associations.collectRefStars(refCat,
1381  self.config.matchCut*lsst.geom.arcseconds,
1382  fluxField,
1383  refCoordinateErr=refCoordErr,
1384  rejectBadFluxes=reject_bad_fluxes)
1385  add_measurement(self.jobjob, 'jointcal.%s_collected_refStars' % name,
1386  associations.refStarListSize())
1387 
1388  associations.prepareFittedStars(self.config.minMeasurements)
1389 
1390  self._check_star_lists_check_star_lists(associations, name)
1391  add_measurement(self.jobjob, 'jointcal.%s_prepared_refStars' % name,
1392  associations.nFittedStarsWithAssociatedRefStar())
1393  add_measurement(self.jobjob, 'jointcal.%s_prepared_fittedStars' % name,
1394  associations.fittedStarListSize())
1395  add_measurement(self.jobjob, 'jointcal.%s_prepared_ccdImages' % name,
1396  associations.nCcdImagesValidForFit())
1397 
1398  load_cat_prof_file = 'jointcal_fit_%s.prof'%name if self.config.detailedProfile else ''
1399  dataName = "{}_{}".format(tract, defaultFilter.physicalLabel)
1400  with pipeBase.cmdLineTask.profile(load_cat_prof_file):
1401  result = fit_function(associations, dataName)
1402  # TODO DM-12446: turn this into a "butler save" somehow.
1403  # Save reference and measurement chi2 contributions for this data
1404  if self.config.writeChi2FilesInitialFinal:
1405  baseName = self._getDebugPath_getDebugPath(f"{name}_final_chi2-{dataName}")
1406  result.fit.saveChi2Contributions(baseName+"{type}")
1407  self.log.info("Wrote chi2 contributions files: %s", baseName)
1408 
1409  return result
1410 
1411  def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radius, filterLabel,
1412  applyColorterms=False, epoch=None):
1413  """Load the necessary reference catalog sources, convert fluxes to
1414  correct units, and apply color term corrections if requested.
1415 
1416  Parameters
1417  ----------
1418  refObjLoader : `lsst.meas.algorithms.LoadReferenceObjectsTask`
1419  The reference catalog loader to use to get the data.
1420  referenceSelector : `lsst.meas.algorithms.ReferenceSourceSelectorTask`
1421  Source selector to apply to loaded reference catalog.
1422  center : `lsst.geom.SpherePoint`
1423  The center around which to load sources.
1424  radius : `lsst.geom.Angle`
1425  The radius around ``center`` to load sources in.
1426  filterLabel : `lsst.afw.image.FilterLabel`
1427  The camera filter to load fluxes for.
1428  applyColorterms : `bool`
1429  Apply colorterm corrections to the refcat for ``filterName``?
1430  epoch : `astropy.time.Time`, optional
1431  Epoch to which to correct refcat proper motion and parallax,
1432  or `None` to not apply such corrections.
1433 
1434  Returns
1435  -------
1436  refCat : `lsst.afw.table.SimpleCatalog`
1437  The loaded reference catalog.
1438  fluxField : `str`
1439  The name of the reference catalog flux field appropriate for ``filterName``.
1440  """
1441  skyCircle = refObjLoader.loadSkyCircle(center,
1442  radius,
1443  filterLabel.bandLabel,
1444  epoch=epoch)
1445 
1446  selected = referenceSelector.run(skyCircle.refCat)
1447  # Need memory contiguity to get reference filters as a vector.
1448  if not selected.sourceCat.isContiguous():
1449  refCat = selected.sourceCat.copy(deep=True)
1450  else:
1451  refCat = selected.sourceCat
1452 
1453  if applyColorterms:
1454  try:
1455  refCatName = refObjLoader.ref_dataset_name # Gen2
1456  except AttributeError:
1457  refCatName = self.config.connections.photometryRefCat # Gen3
1458 
1459  self.log.info("Applying color terms for physical filter=%r reference catalog=%s",
1460  filterLabel.physicalLabel, refCatName)
1461  colorterm = self.config.colorterms.getColorterm(filterLabel.physicalLabel,
1462  refCatName,
1463  doRaise=True)
1464 
1465  refMag, refMagErr = colorterm.getCorrectedMagnitudes(refCat)
1466  refCat[skyCircle.fluxField] = u.Magnitude(refMag, u.ABmag).to_value(u.nJy)
1467  # TODO: I didn't want to use this, but I'll deal with it in DM-16903
1468  refCat[skyCircle.fluxField+'Err'] = fluxErrFromABMagErr(refMagErr, refMag) * 1e9
1469 
1470  return refCat, skyCircle.fluxField
1471 
1472  def _check_star_lists(self, associations, name):
1473  # TODO: these should be len(blah), but we need this properly wrapped first.
1474  if associations.nCcdImagesValidForFit() == 0:
1475  raise RuntimeError('No images in the ccdImageList!')
1476  if associations.fittedStarListSize() == 0:
1477  raise RuntimeError('No stars in the {} fittedStarList!'.format(name))
1478  if associations.refStarListSize() == 0:
1479  raise RuntimeError('No stars in the {} reference star list!'.format(name))
1480 
1481  def _logChi2AndValidate(self, associations, fit, model, chi2Label, writeChi2Name=None):
1482  """Compute chi2, log it, validate the model, and return chi2.
1483 
1484  Parameters
1485  ----------
1486  associations : `lsst.jointcal.Associations`
1487  The star/reference star associations to fit.
1488  fit : `lsst.jointcal.FitterBase`
1489  The fitter to use for minimization.
1490  model : `lsst.jointcal.Model`
1491  The model being fit.
1492  chi2Label : `str`
1493  Label to describe the chi2 (e.g. "Initialized", "Final").
1494  writeChi2Name : `str`, optional
1495  Filename prefix to write the chi2 contributions to.
1496  Do not supply an extension: an appropriate one will be added.
1497 
1498  Returns
1499  -------
1500  chi2: `lsst.jointcal.Chi2Accumulator`
1501  The chi2 object for the current fitter and model.
1502 
1503  Raises
1504  ------
1505  FloatingPointError
1506  Raised if chi2 is infinite or NaN.
1507  ValueError
1508  Raised if the model is not valid.
1509  """
1510  if writeChi2Name is not None:
1511  fullpath = self._getDebugPath_getDebugPath(writeChi2Name)
1512  fit.saveChi2Contributions(fullpath+"{type}")
1513  self.log.info("Wrote chi2 contributions files: %s", fullpath)
1514 
1515  chi2 = fit.computeChi2()
1516  self.log.info("%s %s", chi2Label, chi2)
1517  self._check_stars_check_stars(associations)
1518  if not np.isfinite(chi2.chi2):
1519  raise FloatingPointError(f'{chi2Label} chi2 is invalid: {chi2}')
1520  if not model.validate(associations.getCcdImageList(), chi2.ndof):
1521  raise ValueError("Model is not valid: check log messages for warnings.")
1522  return chi2
1523 
1524  def _fit_photometry(self, associations, dataName=None):
1525  """
1526  Fit the photometric data.
1527 
1528  Parameters
1529  ----------
1530  associations : `lsst.jointcal.Associations`
1531  The star/reference star associations to fit.
1532  dataName : `str`
1533  Name of the data being processed (e.g. "1234_HSC-Y"), for
1534  identifying debugging files.
1535 
1536  Returns
1537  -------
1538  fit_result : `namedtuple`
1539  fit : `lsst.jointcal.PhotometryFit`
1540  The photometric fitter used to perform the fit.
1541  model : `lsst.jointcal.PhotometryModel`
1542  The photometric model that was fit.
1543  """
1544  self.log.info("=== Starting photometric fitting...")
1545 
1546  # TODO: should use pex.config.RegistryField here (see DM-9195)
1547  if self.config.photometryModel == "constrainedFlux":
1548  model = lsst.jointcal.ConstrainedFluxModel(associations.getCcdImageList(),
1549  self.focalPlaneBBoxfocalPlaneBBox,
1550  visitOrder=self.config.photometryVisitOrder,
1551  errorPedestal=self.config.photometryErrorPedestal)
1552  # potentially nonlinear problem, so we may need a line search to converge.
1553  doLineSearch = self.config.allowLineSearch
1554  elif self.config.photometryModel == "constrainedMagnitude":
1555  model = lsst.jointcal.ConstrainedMagnitudeModel(associations.getCcdImageList(),
1556  self.focalPlaneBBoxfocalPlaneBBox,
1557  visitOrder=self.config.photometryVisitOrder,
1558  errorPedestal=self.config.photometryErrorPedestal)
1559  # potentially nonlinear problem, so we may need a line search to converge.
1560  doLineSearch = self.config.allowLineSearch
1561  elif self.config.photometryModel == "simpleFlux":
1562  model = lsst.jointcal.SimpleFluxModel(associations.getCcdImageList(),
1563  errorPedestal=self.config.photometryErrorPedestal)
1564  doLineSearch = False # purely linear in model parameters, so no line search needed
1565  elif self.config.photometryModel == "simpleMagnitude":
1566  model = lsst.jointcal.SimpleMagnitudeModel(associations.getCcdImageList(),
1567  errorPedestal=self.config.photometryErrorPedestal)
1568  doLineSearch = False # purely linear in model parameters, so no line search needed
1569 
1570  fit = lsst.jointcal.PhotometryFit(associations, model)
1571  # TODO DM-12446: turn this into a "butler save" somehow.
1572  # Save reference and measurement chi2 contributions for this data
1573  if self.config.writeChi2FilesInitialFinal:
1574  baseName = f"photometry_initial_chi2-{dataName}"
1575  else:
1576  baseName = None
1577  if self.config.writeInitialModel:
1578  fullpath = self._getDebugPath_getDebugPath(f"initial_photometry_model-{dataName}.txt")
1579  writeModel(model, fullpath, self.log)
1580  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialized", writeChi2Name=baseName)
1581 
1582  def getChi2Name(whatToFit):
1583  if self.config.writeChi2FilesOuterLoop:
1584  return f"photometry_init-%s_chi2-{dataName}" % whatToFit
1585  else:
1586  return None
1587 
1588  # The constrained model needs the visit transform fit first; the chip
1589  # transform is initialized from the singleFrame PhotoCalib, so it's close.
1590  if self.config.writeInitMatrix:
1591  dumpMatrixFile = self._getDebugPath_getDebugPath(f"photometry_preinit-{dataName}")
1592  else:
1593  dumpMatrixFile = ""
1594  if self.config.photometryModel.startswith("constrained"):
1595  # no line search: should be purely (or nearly) linear,
1596  # and we want a large step size to initialize with.
1597  fit.minimize("ModelVisit", dumpMatrixFile=dumpMatrixFile)
1598  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialize ModelVisit",
1599  writeChi2Name=getChi2Name("ModelVisit"))
1600  dumpMatrixFile = "" # so we don't redo the output on the next step
1601 
1602  fit.minimize("Model", doLineSearch=doLineSearch, dumpMatrixFile=dumpMatrixFile)
1603  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialize Model",
1604  writeChi2Name=getChi2Name("Model"))
1605 
1606  fit.minimize("Fluxes") # no line search: always purely linear.
1607  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialize Fluxes",
1608  writeChi2Name=getChi2Name("Fluxes"))
1609 
1610  fit.minimize("Model Fluxes", doLineSearch=doLineSearch)
1611  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialize ModelFluxes",
1612  writeChi2Name=getChi2Name("ModelFluxes"))
1613 
1614  model.freezeErrorTransform()
1615  self.log.debug("Photometry error scales are frozen.")
1616 
1617  chi2 = self._iterate_fit_iterate_fit(associations,
1618  fit,
1619  self.config.maxPhotometrySteps,
1620  "photometry",
1621  "Model Fluxes",
1622  doRankUpdate=self.config.photometryDoRankUpdate,
1623  doLineSearch=doLineSearch,
1624  dataName=dataName)
1625 
1626  add_measurement(self.jobjob, 'jointcal.photometry_final_chi2', chi2.chi2)
1627  add_measurement(self.jobjob, 'jointcal.photometry_final_ndof', chi2.ndof)
1628  return Photometry(fit, model)
1629 
1630  def _fit_astrometry(self, associations, dataName=None):
1631  """
1632  Fit the astrometric data.
1633 
1634  Parameters
1635  ----------
1636  associations : `lsst.jointcal.Associations`
1637  The star/reference star associations to fit.
1638  dataName : `str`
1639  Name of the data being processed (e.g. "1234_HSC-Y"), for
1640  identifying debugging files.
1641 
1642  Returns
1643  -------
1644  fit_result : `namedtuple`
1645  fit : `lsst.jointcal.AstrometryFit`
1646  The astrometric fitter used to perform the fit.
1647  model : `lsst.jointcal.AstrometryModel`
1648  The astrometric model that was fit.
1649  sky_to_tan_projection : `lsst.jointcal.ProjectionHandler`
1650  The model for the sky to tangent plane projection that was used in the fit.
1651  """
1652 
1653  self.log.info("=== Starting astrometric fitting...")
1654 
1655  associations.deprojectFittedStars()
1656 
1657  # NOTE: need to return sky_to_tan_projection so that it doesn't get garbage collected.
1658  # TODO: could we package sky_to_tan_projection and model together so we don't have to manage
1659  # them so carefully?
1660  sky_to_tan_projection = lsst.jointcal.OneTPPerVisitHandler(associations.getCcdImageList())
1661 
1662  if self.config.astrometryModel == "constrained":
1663  model = lsst.jointcal.ConstrainedAstrometryModel(associations.getCcdImageList(),
1664  sky_to_tan_projection,
1665  chipOrder=self.config.astrometryChipOrder,
1666  visitOrder=self.config.astrometryVisitOrder)
1667  elif self.config.astrometryModel == "simple":
1668  model = lsst.jointcal.SimpleAstrometryModel(associations.getCcdImageList(),
1669  sky_to_tan_projection,
1670  self.config.useInputWcs,
1671  nNotFit=0,
1672  order=self.config.astrometrySimpleOrder)
1673 
1674  fit = lsst.jointcal.AstrometryFit(associations, model, self.config.positionErrorPedestal)
1675  # TODO DM-12446: turn this into a "butler save" somehow.
1676  # Save reference and measurement chi2 contributions for this data
1677  if self.config.writeChi2FilesInitialFinal:
1678  baseName = f"astrometry_initial_chi2-{dataName}"
1679  else:
1680  baseName = None
1681  if self.config.writeInitialModel:
1682  fullpath = self._getDebugPath_getDebugPath(f"initial_astrometry_model-{dataName}.txt")
1683  writeModel(model, fullpath, self.log)
1684  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initial", writeChi2Name=baseName)
1685 
1686  def getChi2Name(whatToFit):
1687  if self.config.writeChi2FilesOuterLoop:
1688  return f"astrometry_init-%s_chi2-{dataName}" % whatToFit
1689  else:
1690  return None
1691 
1692  if self.config.writeInitMatrix:
1693  dumpMatrixFile = self._getDebugPath_getDebugPath(f"astrometry_preinit-{dataName}")
1694  else:
1695  dumpMatrixFile = ""
1696  # The constrained model needs the visit transform fit first; the chip
1697  # transform is initialized from the detector's cameraGeom, so it's close.
1698  if self.config.astrometryModel == "constrained":
1699  fit.minimize("DistortionsVisit", dumpMatrixFile=dumpMatrixFile)
1700  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialize DistortionsVisit",
1701  writeChi2Name=getChi2Name("DistortionsVisit"))
1702  dumpMatrixFile = "" # so we don't redo the output on the next step
1703 
1704  fit.minimize("Distortions", dumpMatrixFile=dumpMatrixFile)
1705  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialize Distortions",
1706  writeChi2Name=getChi2Name("Distortions"))
1707 
1708  fit.minimize("Positions")
1709  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialize Positions",
1710  writeChi2Name=getChi2Name("Positions"))
1711 
1712  fit.minimize("Distortions Positions")
1713  self._logChi2AndValidate_logChi2AndValidate(associations, fit, model, "Initialize DistortionsPositions",
1714  writeChi2Name=getChi2Name("DistortionsPositions"))
1715 
1716  chi2 = self._iterate_fit_iterate_fit(associations,
1717  fit,
1718  self.config.maxAstrometrySteps,
1719  "astrometry",
1720  "Distortions Positions",
1721  sigmaRelativeTolerance=self.config.astrometryOutlierRelativeTolerance,
1722  doRankUpdate=self.config.astrometryDoRankUpdate,
1723  dataName=dataName)
1724 
1725  add_measurement(self.jobjob, 'jointcal.astrometry_final_chi2', chi2.chi2)
1726  add_measurement(self.jobjob, 'jointcal.astrometry_final_ndof', chi2.ndof)
1727 
1728  return Astrometry(fit, model, sky_to_tan_projection)
1729 
1730  def _check_stars(self, associations):
1731  """Count measured and reference stars per ccd and warn/log them."""
1732  for ccdImage in associations.getCcdImageList():
1733  nMeasuredStars, nRefStars = ccdImage.countStars()
1734  self.log.debug("ccdImage %s has %s measured and %s reference stars",
1735  ccdImage.getName(), nMeasuredStars, nRefStars)
1736  if nMeasuredStars < self.config.minMeasuredStarsPerCcd:
1737  self.log.warning("ccdImage %s has only %s measuredStars (desired %s)",
1738  ccdImage.getName(), nMeasuredStars, self.config.minMeasuredStarsPerCcd)
1739  if nRefStars < self.config.minRefStarsPerCcd:
1740  self.log.warning("ccdImage %s has only %s RefStars (desired %s)",
1741  ccdImage.getName(), nRefStars, self.config.minRefStarsPerCcd)
1742 
1743  def _iterate_fit(self, associations, fitter, max_steps, name, whatToFit,
1744  dataName="",
1745  sigmaRelativeTolerance=0,
1746  doRankUpdate=True,
1747  doLineSearch=False):
1748  """Run fitter.minimize up to max_steps times, returning the final chi2.
1749 
1750  Parameters
1751  ----------
1752  associations : `lsst.jointcal.Associations`
1753  The star/reference star associations to fit.
1754  fitter : `lsst.jointcal.FitterBase`
1755  The fitter to use for minimization.
1756  max_steps : `int`
1757  Maximum number of steps to run outlier rejection before declaring
1758  convergence failure.
1759  name : {'photometry' or 'astrometry'}
1760  What type of data are we fitting (for logs and debugging files).
1761  whatToFit : `str`
1762  Passed to ``fitter.minimize()`` to define the parameters to fit.
1763  dataName : `str`, optional
1764  Descriptive name for this dataset (e.g. tract and filter),
1765  for debugging.
1766  sigmaRelativeTolerance : `float`, optional
1767  Convergence tolerance for the fractional change in the chi2 cut
1768  level for determining outliers. If set to zero, iterations will
1769  continue until there are no outliers.
1770  doRankUpdate : `bool`, optional
1771  Do an Eigen rank update during minimization, or recompute the full
1772  matrix and gradient?
1773  doLineSearch : `bool`, optional
1774  Do a line search for the optimum step during minimization?
1775 
1776  Returns
1777  -------
1778  chi2: `lsst.jointcal.Chi2Statistic`
1779  The final chi2 after the fit converges, or is forced to end.
1780 
1781  Raises
1782  ------
1783  FloatingPointError
1784  Raised if the fitter fails with a non-finite value.
1785  RuntimeError
1786  Raised if the fitter fails for some other reason;
1787  log messages will provide further details.
1788  """
1789  if self.config.writeInitMatrix:
1790  dumpMatrixFile = self._getDebugPath_getDebugPath(f"{name}_postinit-{dataName}")
1791  else:
1792  dumpMatrixFile = ""
1793  oldChi2 = lsst.jointcal.Chi2Statistic()
1794  oldChi2.chi2 = float("inf")
1795  for i in range(max_steps):
1796  if self.config.writeChi2FilesOuterLoop:
1797  writeChi2Name = f"{name}_iterate_{i}_chi2-{dataName}"
1798  else:
1799  writeChi2Name = None
1800  result = fitter.minimize(whatToFit,
1801  self.config.outlierRejectSigma,
1802  sigmaRelativeTolerance=sigmaRelativeTolerance,
1803  doRankUpdate=doRankUpdate,
1804  doLineSearch=doLineSearch,
1805  dumpMatrixFile=dumpMatrixFile)
1806  dumpMatrixFile = "" # clear it so we don't write the matrix again.
1807  chi2 = self._logChi2AndValidate_logChi2AndValidate(associations, fitter, fitter.getModel(),
1808  f"Fit iteration {i}", writeChi2Name=writeChi2Name)
1809 
1810  if result == MinimizeResult.Converged:
1811  if doRankUpdate:
1812  self.log.debug("fit has converged - no more outliers - redo minimization "
1813  "one more time in case we have lost accuracy in rank update.")
1814  # Redo minimization one more time in case we have lost accuracy in rank update
1815  result = fitter.minimize(whatToFit, self.config.outlierRejectSigma,
1816  sigmaRelativeTolerance=sigmaRelativeTolerance)
1817  chi2 = self._logChi2AndValidate_logChi2AndValidate(associations, fitter, fitter.getModel(), "Fit completed")
1818 
1819  # log a message for a large final chi2, TODO: DM-15247 for something better
1820  if chi2.chi2/chi2.ndof >= 4.0:
1821  self.log.error("Potentially bad fit: High chi-squared/ndof.")
1822 
1823  break
1824  elif result == MinimizeResult.Chi2Increased:
1825  self.log.warning("Still some outliers remaining but chi2 increased - retry")
1826  # Check whether the increase was large enough to cause trouble.
1827  chi2Ratio = chi2.chi2 / oldChi2.chi2
1828  if chi2Ratio > 1.5:
1829  self.log.warning('Significant chi2 increase by a factor of %.4g / %.4g = %.4g',
1830  chi2.chi2, oldChi2.chi2, chi2Ratio)
1831  # Based on a variety of HSC jointcal logs (see DM-25779), it
1832  # appears that chi2 increases more than a factor of ~2 always
1833  # result in the fit diverging rapidly and ending at chi2 > 1e10.
1834  # Using 10 as the "failure" threshold gives some room between
1835  # leaving a warning and bailing early.
1836  if chi2Ratio > 10:
1837  msg = ("Large chi2 increase between steps: fit likely cannot converge."
1838  " Try setting one or more of the `writeChi2*` config fields and looking"
1839  " at how individual star chi2-values evolve during the fit.")
1840  raise RuntimeError(msg)
1841  oldChi2 = chi2
1842  elif result == MinimizeResult.NonFinite:
1843  filename = self._getDebugPath_getDebugPath("{}_failure-nonfinite_chi2-{}.csv".format(name, dataName))
1844  # TODO DM-12446: turn this into a "butler save" somehow.
1845  fitter.saveChi2Contributions(filename+"{type}")
1846  msg = "Nonfinite value in chi2 minimization, cannot complete fit. Dumped star tables to: {}"
1847  raise FloatingPointError(msg.format(filename))
1848  elif result == MinimizeResult.Failed:
1849  raise RuntimeError("Chi2 minimization failure, cannot complete fit.")
1850  else:
1851  raise RuntimeError("Unxepected return code from minimize().")
1852  else:
1853  self.log.error("%s failed to converge after %d steps"%(name, max_steps))
1854 
1855  return chi2
1856 
1857  def _make_output(self, ccdImageList, model, func):
1858  """Return the internal jointcal models converted to the afw
1859  structures that will be saved to disk.
1860 
1861  Parameters
1862  ----------
1863  ccdImageList : `lsst.jointcal.CcdImageList`
1864  The list of CcdImages to get the output for.
1865  model : `lsst.jointcal.AstrometryModel` or `lsst.jointcal.PhotometryModel`
1866  The internal jointcal model to convert for each `lsst.jointcal.CcdImage`.
1867  func : `str`
1868  The name of the function to call on ``model`` to get the converted
1869  structure. Must accept an `lsst.jointcal.CcdImage`.
1870 
1871  Returns
1872  -------
1873  output : `dict` [`tuple`, `lsst.jointcal.AstrometryModel`] or
1874  `dict` [`tuple`, `lsst.jointcal.PhotometryModel`]
1875  The data to be saved, keyed on (visit, detector).
1876  """
1877  output = {}
1878  for ccdImage in ccdImageList:
1879  ccd = ccdImage.ccdId
1880  visit = ccdImage.visit
1881  self.log.debug("%s for visit: %d, ccd: %d", func, visit, ccd)
1882  output[(visit, ccd)] = getattr(model, func)(ccdImage)
1883  return output
1884 
1885  def _write_astrometry_results(self, associations, model, visit_ccd_to_dataRef):
1886  """
1887  Write the fitted astrometric results to a new 'jointcal_wcs' dataRef.
1888 
1889  Parameters
1890  ----------
1891  associations : `lsst.jointcal.Associations`
1892  The star/reference star associations to fit.
1893  model : `lsst.jointcal.AstrometryModel`
1894  The astrometric model that was fit.
1895  visit_ccd_to_dataRef : `dict` of Key: `lsst.daf.persistence.ButlerDataRef`
1896  Dict of ccdImage identifiers to dataRefs that were fit.
1897  """
1898  ccdImageList = associations.getCcdImageList()
1899  output = self._make_output_make_output(ccdImageList, model, "makeSkyWcs")
1900  for key, skyWcs in output.items():
1901  dataRef = visit_ccd_to_dataRef[key]
1902  try:
1903  dataRef.put(skyWcs, 'jointcal_wcs')
1904  except pexExceptions.Exception as e:
1905  self.log.fatal('Failed to write updated Wcs: %s', str(e))
1906  raise e
1907 
1908  def _write_photometry_results(self, associations, model, visit_ccd_to_dataRef):
1909  """
1910  Write the fitted photometric results to a new 'jointcal_photoCalib' dataRef.
1911 
1912  Parameters
1913  ----------
1914  associations : `lsst.jointcal.Associations`
1915  The star/reference star associations to fit.
1916  model : `lsst.jointcal.PhotometryModel`
1917  The photoometric model that was fit.
1918  visit_ccd_to_dataRef : `dict` of Key: `lsst.daf.persistence.ButlerDataRef`
1919  Dict of ccdImage identifiers to dataRefs that were fit.
1920  """
1921 
1922  ccdImageList = associations.getCcdImageList()
1923  output = self._make_output_make_output(ccdImageList, model, "toPhotoCalib")
1924  for key, photoCalib in output.items():
1925  dataRef = visit_ccd_to_dataRef[key]
1926  try:
1927  dataRef.put(photoCalib, 'jointcal_photoCalib')
1928  except pexExceptions.Exception as e:
1929  self.log.fatal('Failed to write updated PhotoCalib: %s', str(e))
1930  raise e
table::Key< int > type
Definition: Detector.cc:163
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:382
static Schema makeMinimalSchema()
Return a minimal schema for Source tables and records.
Definition: Source.h:258
Class for storing ordered metadata with comments.
Definition: PropertyList.h:68
A class representing an angle.
Definition: Angle.h:127
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.
Definition: Associations.h:54
Class that handles the astrometric least squares problem.
Definition: AstrometryFit.h:78
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.
Definition: PhotometryFit.h:45
A model where there is one independent transform per CcdImage.
def getTargetList(parsedCmd, **kwargs)
Definition: jointcal.py:78
def _load_data(self, inputSourceTableVisit, inputVisitSummary, associations, jointcalControl, camera)
Definition: jointcal.py:925
def _compute_proper_motion_epoch(self, ccdImageList)
Definition: jointcal.py:1308
def _get_refcat_coordinate_error_override(self, refCat, name)
Definition: jointcal.py:1261
def _getDebugPath(self, filename)
Definition: jointcal.py:1159
def _check_star_lists(self, associations, name)
Definition: jointcal.py:1472
def _make_one_input_data(self, visitRecord, catalog, detectorDict)
Definition: jointcal.py:989
def _iterate_fit(self, associations, fitter, max_steps, name, whatToFit, dataName="", sigmaRelativeTolerance=0, doRankUpdate=True, doLineSearch=False)
Definition: jointcal.py:1747
def _write_astrometry_results(self, associations, model, visit_ccd_to_dataRef)
Definition: jointcal.py:1885
def _check_stars(self, associations)
Definition: jointcal.py:1730
def _prep_sky(self, associations, filters)
Definition: jointcal.py:1164
def _readDataId(self, butler, dataId)
Definition: jointcal.py:1109
def runQuantum(self, butlerQC, inputRefs, outputRefs)
Definition: jointcal.py:693
def _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:1329
def _put_metrics(self, butlerQC, job, outputRefs)
Definition: jointcal.py:722
def _build_ccdImage(self, data, associations, jointcalControl)
Definition: jointcal.py:1064
def runDataRef(self, dataRefs)
Definition: jointcal.py:1187
def __init__(self, butler=None, initInputs=None, **kwargs)
Definition: jointcal.py:672
def _get_sourceTable_visit_columns(self, inColumns)
Definition: jointcal.py:1003
def run(self, inputSourceTableVisit, inputVisitSummary, inputCamera, tract=None)
Definition: jointcal.py:797
def _fit_photometry(self, associations, dataName=None)
Definition: jointcal.py:1524
def _write_photometry_results(self, associations, model, visit_ccd_to_dataRef)
Definition: jointcal.py:1908
def _logChi2AndValidate(self, associations, fit, model, chi2Label, writeChi2Name=None)
Definition: jointcal.py:1481
def _extract_detector_catalog_from_visit_catalog(self, table, visitCatalog, detectorId, detectorColumn, ixxColumns)
Definition: jointcal.py:875
def _make_output(self, ccdImageList, model, func)
Definition: jointcal.py:1857
def _put_output(self, butlerQC, outputs, outputRefs, camera, setter)
Definition: jointcal.py:738
def _load_reference_catalog(self, refObjLoader, referenceSelector, center, radius, filterLabel, applyColorterms=False, epoch=None)
Definition: jointcal.py:1412
def _fit_astrometry(self, associations, dataName=None)
Definition: jointcal.py:1630
def loadData(self, dataRefs, associations, jointcalControl)
Definition: jointcal.py:1135
Provides consistent interface for LSST exceptions.
Definition: Exception.h:107
std::shared_ptr< FrameSet > append(FrameSet const &first, FrameSet const &second)
Construct a FrameSet that performs two transformations in series.
Definition: functional.cc:33
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
double fluxErrFromABMagErr(double magErr, double mag) noexcept
Compute flux error in Janskys from AB magnitude error and AB magnitude.
Definition: Calib.h:60
def lookupStaticCalibrations(datasetType, registry, quantumDataId, collections)
Definition: jointcal.py:142
def add_measurement(job, name, value)
Definition: jointcal.py:60
def lookupVisitRefCats(datasetType, registry, quantumDataId, collections)
Definition: jointcal.py:182
def writeModel(model, filename, log)
Definition: jointcal.py:626
Definition: Log.h:717
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
def format(config, name=None, writeSourceLine=True, prefix="", verbose=False)
Definition: history.py:174