How to use WCSs in DP1 and commissioning processing

Support Data Products
, ,
#1

We’ve seen a lot of confusion with our WCS objects recently, in both the DP1 release and with commissioning scientists encountering them in internal processing. It’s become clear that some aspects of the lsst.afw.geom.SkyWcs interface are (to be generous) super confusing, and since it’ll take it a while to fix that directly (see RFC-1106 if you’re interested in the details of that), I wanted to document here what the WCSs on various data products actually are, and which operations behave more-or-less as one would reasonably expect (and which do not).

  • The WCS objects attached to our raw datasets (i.e. what you get from butler.get("raw", ...).getWcs()) are formed by composing a static estimate of the camera geometry with an on-sky shift and rotation from the boresight angle and rotation angle in the headers. The biggest problem with these for LSSTCam and ComCam is that the boresight pointing can be pretty inaccurate. In addition, for both instruments, the static camera geometry is something we guessed at before the instruments went on sky, not something we fit afterwards (that will change in the near future for LSSTCam), and we know it’s not quite right. There may also be a FITS WCS in the raw headers coming from the instrument, but DM does not ever use it, and I can’t speak to its quality at all.

  • The WCS objects attached to our preliminary_visit_image and preliminary_visit_summary datasets (not included in DP1) are from a simple affine transform fit that maps the raw WCS positions to our Gaia-based reference catalog. This fully addresses the boresight pointing problem and mitigates the issues with the static camera geometry, at least on larger scales and when our detection of calibration stars worked well (it notably does not work well in crowded fields right now). This WCS is not representable as FITS, and at present if you think you see a FITS WCS version of it somewhere, it is quite likely to be garbage (more on this later).

  • The WCS objects attached to our visit_image and visit_summary datasets are from a moderately sophisticated multi-epoch fit using GBDES, the same code used to produce some excellent DECam astrometry (but, at present, without the pixel-level and Gaussian process terms that made that work especially excellent). To our knowledge these WCSs are quite good, and should get ever better in the future. However, they are also not representable as FITS, and it is all too easy to get a FITS representation that is garbage. For DP1 only, we have attached an explicit TAN-SIP approximation to the visit_image (only!) that should be good to about 0.1 pixel over the full image, with the idea that this should be good enough for visualization. This is what’s stored in the visit_image FITS headers. But outside DP1, or if you get a DP1 FITS representation via something other than the FITS file headers or by calling the wcs.getFitsApproximation() method, it is likely to be garbage.

  • The WCS objects attached to our coadds come from a predefined tiling of the sky (not some kind of fitting), and are fully FITS representable (they’re all strictly gnomonic, i.e. TAN). It is of course possible for there to be problems in the WCSs of the visits that went into the coadds, which then looks like an astrometric problem with the coadds themselves (this is known to be the case at least sometimes in crowded fields). For coadd WCSs, the big gotcha is that they correspond to a full-tract coordinate system, which is offset from the image array indexes for any particular patch by image.getXY0() or image.getBBox().getMin(). The same would be true for the WCS attached to any cutout of a [preliminary_]visit_image). We conceptually think of this as the first pixel in the image having an index that is not in general (0, 0).

When you have a lsst.afw.geom.SkyWcs object (what exposure.wcs is), generally the only things you should do with it are:

  • call pixelToSky[Array] to transform points from pixel coordinates (which may not start at (0, 0)!) to celestial coordinates;
  • call skyToPixel[Array] to do the reverse;
  • call linearizePixelToSky and linearizeSkyToPixel to get local affine approximations to the WCS at a point.
  • call getPixelScale(point) (with the position to evaluate at!) to get the local pixel scale.

For coadd WCSs only, wcs.getFitsMetadata().toDict() is also safe to call, as a way to get the header keys for the corresponding FITS WCS. And for DP1 visit_image and visit_summary WCS objects, wcs.getFitsApproximation().getFitsMetadata().toDict() is safe (but an approximation!). I’m afraid in other cases there just isn’t a FITS approximation available; we’re working to rectify this, but in the meantime you’re stuck either using our WCS objects directly or fitting your own FITS approximation to point-pairs you’ve asked ours to generate.

Finally, please do not:

  • call getPixelOrigin(). This might be CRPIX for a FITS WCS, or it might be garbage;
  • call getSkyOrigin(). This might be CRVAL for a FITS WCS, or it might be garbage;
  • call getCdMatrix(). This might be the CD matrix for a FITS WCS, or it might be garbage;
  • call getPixelScale() (without arguments). This is evaluated at getPixelOrigin(), and hence may be garbage. If you pass a point in explicitly, this is fine.
  • call getTanWcs(). This gives you a FITS-compatible TAN WCS, but it might be at a reference point far from the actual image you’re looking at and it might have little in common with the WCS you started with. It’s likely to look semi-okay, in fact, for most of the above cases, which actually just makes it more dangerous.

As you can see in the RFC linked above, I’ve proposed removing all of these methods.

1 Like