Jointcal is a package for performing simultaneous astrometric and photometric fitting of multiple visits of the same field to determine the optical distortions and photometric calibration of the system as a whole. I have run it on a variety of test data, but it’s time to start poking at it in other ways, hence this beta test announcement.
Requirements
You can install jointcal in the usual manner (e.g., rebuild jointcal in lsstsw).
optional dependency on testdata_jointcal (a little over 1 GB of pre-processed catalogs from multiple instruments) in order to run the included unittests.
a butler repository processed with stack v13.0 or higher (jointcal requires processed catalogs with a fully-populated VisitInfo).
Usage
Jointcal can be executed as a CmdlineTask in the usual manner, specifying an input butler repository, output directory, and data ids. See examples/jointcal_validation_data_hsc.sh for an example of how to run jointcal as a CmdlineTask on the validation_data_hsc dataset (Note: that product is ~400GB git lfs; you’ll have to reprocess it to v13.0 or wait for it to be updated to that level, in order to use that script).
Configurables of interest include astrometryRefObjLoader (defaults to Astrometry.net, but can be retargetted to LoadIndexedReferenceObjectsTask) and doPhotometry (given caveat 2. below, you may want to set this to “False”).
The fitted model is written to a butler “wcs” Exposure object, containing the fitted Wcs and Calib, but no pixels or other metata. One can use lsst.afw.table.utils.updateSourceCoords(wcs.getWcs(), cat) to update in-place the Coords in a Catalog, given that new wcs. Jointcal also includes bin/plot_jointcal_results.py, to help evaluate the results of the fit.
Caveats:
The default astrometric fitter does not include the full-focal plane term (see DM-9506 for the plan to make that configurable). Even with that, the astrometric fit should be noticeably better than single frame processing, and will improve with more visits of a field being available.
The photometric fitter will do no better than single frame processing (see DM-9184 for progress on the next step).
The fitting step is very fast, but a fair bit of the time is taken up reading the calexps (which should improve drastically once DM-9153 is completed).
Once points 1. and 2. are dealt with, jointcal will be RFC’d to become part of lsst_distrib. For the HSC folks out there, that’s also the point at which I’d say we could probably replace meas_mosaic with jointcal (possibly after some head-to-head comparisons).
Thanks to @boutigny and @PierreAstier for their initial work on meas_simastrom, on which jointcal is based.
I’m working on the warping code, and would love to include a hook to use the results of jointcal. Could you please outline the code required to get a corrected calexp (i.e., with jointcal calibrations applied instead of the ones from single frame measurement)?
If you want to update the catalogs in a list of dataRefs that have been processed with jointcal, this should do the trick:
wcss = [ref.get('wcs') for ref in dataRefs]
catalogs = [ref.get('src') for ref in dataRefs]
for wcs, cat in zip(wcss, catalogs):
lsst.afw.table.utils.updateSourceCoords(wcs.getWcs(), cat)
Images? I’m confused. Do you mean replacing the wcs in the calexp with the wcs fitted by jointcal? This should work:
wcss = [ref.get('wcs') for ref in dataRefs]
calexps = [ref.get('calexp') for ref in dataRefs]
for wcs, calexp in zip(wcss, calexps):
calexp.setWcs(wcs.getWcs())
Obviously, you’d have to butler put after this, if you’re doing this in other code (same with the snippet above).
The wcs is one element that you’re correcting. The other is the flux (currently through a single zero point per CCD, but you said that’s going to get more complicated). I suggest you should have a single function to retrieve a corrected calexp (e.g., as in meas_mosaic).
Good point. I’m recommending ignoring the fitted Calib right now, because the current model hasn’t produced results measureably different from singleFrame in any of my tests (see caveat 2).
I’m torn about adding a method to update a calexp in-place: the butler redesign is supposed to handle that “automatically”. Perhaps @natepease or @ktl can provide some insight on the question of composite datasets that have components updated by subsequent processing (that live in a different dataset).
If you have a composite calib dataset that stored the wcs separately, you should be able to put an updated wcs dataset and then get the calib, and the wcs that the composite gets assembled with would be the newer of the wcs datasets. (Does that answer the insight you were looking for?)
Regardless of what the current model is (single CCD or polynomial) or where it happens (the butler or science pipelines code), the calexp is going to need to be updated, and the code that knows how to do it is in jointcal, because it’s jointcal that provided the persisted model in the first place. That means that the code to do the update must live in jointcal somewhere.
That’s true in the meas_mosaic model, where at least some of the outputs are rather weird things we don’t expect any other code to understand. We’d like to eventually get the jointcal results into regular afw objects we expect all code to understand (e.g. Calib). As for the Wcs, I believe the jointcal outputs can already be used as a drop-in replacement for the meas_mosaic ones. If the code to apply those wasn’t tied up with the code to apply the photometric calibrations, I’d advocating moving it out of meas_mosaic.