Incorrect WCS returned by getWcs().getFitsMetadata() for cutouts?

Support
,
#1

If I create a cutout (either from a coadd by following the instructions in https://github.com/lsst/tutorial-notebooks/blob/main/DP0.2/03c_Big_deepCoadd_Cutout.ipynb or from single exposures by following tutorial-notebooks/DP0.2/13a_Image_Cutout_SciDemo.ipynb at main · lsst/tutorial-notebooks · GitHub), I find that the WCS returned by running, for example, wcs_g = WCS(cutout_image_g.getWcs().getFitsMetadata()) has different (and incorrect) CRPIX and NAXIS values compared to when I instead do

filename = os.path.join(tempdir,'cutout_g.fits')
cutout_image_g.writeFits(filename)
f=fits.open(filename)
wcs_g = WCS(f[1].header)

which successfully produces the correct WCS. I was wondering if this is expected for cutouts, and if there’s a neater way to extract the correct WCS without saving out to a temp file?

Many thanks!

#2

The issue is that our WCS objects do not hold the offset of the origin of the image they are attached to; we prefer to put that offset into the image itself only. This saves us from having to create slightly-modified WCS objects whenever we make a subimage, and it allows a lot more code to not care about whether an image has a WCS at all. And putting the offset in both the image and the WCS would open the door to those getting out of sync.

But FITS took the opposite approach: it has no concept of what the origin of an image is (it always assumes it’s (1, 1) - which is notably different from the (0, 0) most Python codes assume), and so we have to move the offset from our image objects into the WCS representation when we write to FITS.

You can also do that yourself, though it’s a little verbose:

from lsst.geom import Point2I

shifted_wcs = cutout_image.getWcs().copyAtShiftedPixelOrigin(Point2I(0, 0) - cutout_image.getXY0())
print(shifted_wcs.getFitsMetadata())
#3

That makes sense, thank you! A follow-up question, copyAtShiftedPixelOrigin() seems to expect two arguments: (self: lsst.afw.geom.SkyWcs, shift: lsst.geom.Extent2D). I assume the first argument is just the cutout wcs, but for the second argument, is there a way to convert the Point2I(0, 0) - cutout_image.getXY0() argument to an Extent2D argument? Many thanks!

#4

Ah, oops, sorry about the bad example. Hopefully I won’t get it wrong twice:

from lsst.geom import Extent2D, Point2I

offset = Extent2D(Point2I(0, 0) - cutout_image.getXY0())
shifted_wcs = cutout_image.getWcs().copyAtShiftedPixelOrigin(cutout.getWcs(), offset)
#5

Thank you very much! This worked well with Extent2D (and it turns out that’s the only argument needed for copyAtShiftedPixelOrigin. Just for reference for others the full solution in the end was:

import lsst.geom as geom
from astropy.wcs import WCS
offset_g = geom.Extent2D(geom.Point2I(0, 0) - cutout_image_g.getXY0())
shifted_wcs_g = cutout_image_g.getWcs().copyAtShiftedPixelOrigin(offset)
wcs_g = WCS(shifted_wcs_g.getFitsMetadata())