Rubin image WCS sky coordinates to pixel position

Support Data Products
#1

Hi, I have been playing with querying the images associated with various diaSources and I have noticed some offsets in the pixel coordinates for a given sky coordinate depending on which WCS I use.

The ExposureF images have a wcs attribute which I understand to be the most accurate representation of the image WCS. I have found examples of how this WCS can be approximated in FITS file format (e.g. by converting the ExposureF WCS or when saving a cutout to a FITS file). It is stated that the FITS approximation is of course not exact but that it should be accurate to within 0.05 pixels.

I might be doing something wrong but I’ve tested a hnadful of visits and I’m finding positional offset of ~5 pixels when I try the two different WCS methods, which is much bigger than I was expecting. Any comments or advice? Here is a demo notebook: https://github.com/jrob93/lsst_nbs/blob/main/dp1/diaSource_WCS_analysis.ipynb

#2

Hi @jrob93 thanks for your question. In your notebook you used skycoord_to_pixel from astropy.wcs.utils. There could be small definition differences between Astropy and the visit_images. Perhaps you could try the following scripts. Let’s wait and see if DM experts can comment on this.

import lsst.geom as geom

x, y = wcs.skyToPixel(geom.SpherePoint(ra, dec, geom.degrees))

from lsst.daf.butler import Butler

butler = Butler("dp1", collections="LSSTComCam/DP1")
query = "band.name = :band AND visit_detector_region.region OVERLAPS POINT(:ra, :dec)"
bind = {'band': band, 'ra': ra, 'dec': dec}
visit_image_refs = butler.query_datasets("visit_image", where=query,
                                         bind=bind, order_by=["visit.timespan.begin"])
visit_image = butler.get(visit_image_refs[0])

fig = plt.figure()
plt.subplot(projection=WCS(visit_image.getWcs().getFitsMetadata()))
1 Like
#3

That way of extracting a FITS WCS is indeed unsafe for visit_image or difference_image; instead of

wcs_fits = WCS(exp.getWcs().getFitsMetadata())

you need to do

wcs_fits = WCS(exp.getWcs().getFitsApproximation().getFitsMetadata())

The former is safe only for coadd WCSs, which are natively FITS-compatible.

#4

Thanks both for the replies. @jbosch I’ve tried your suggestion for the safe wcs_fits but it seems to result in the exact same WCS and pixel offsets (at least for the handful of example visits I’m looking at). E.g. for visit = 2024112700171:

Screenshot 2025-12-01 at 20.13.24
Screenshot 2025-12-01 at 20.13.241232Ă—942 123 KB

#5

Sorry, I didn’t remember how things worked in detail and didn’t take the time to check; getFitsMetadata will call getFitsApproximation internally in this case and so what you were doing is safe, and my recommendation was not at all helpful.

I’ve managed to reproduce your result in the attached notebook (this was running on an internal copy of the data, so you might need to adjust the Butler initialization slightly to run it):

DP1-SIP-debugging.ipynb (9.7 KB)

Our SIP approximation (as evaluated by our code) differs from our “true” WCS in the way we’ve previously reported. But an Astropy WCS initialized with the same SIP coefficients (and I’ve checked that they are the same in the Astropy WCS object) transforms points differently, especially in the sky-to-pixel direction, where the differences are more than a pixel.

I’m pretty confident in our TAN-SIP evaluation logic - we delegate basically all of this to the StarLink AST library - but a bug like this in Astropy seems pretty unlikely, too. This will take some digging through both the Astropy code and ours to figure out what’s going on. I’ll also see if I can find any other independent implementations to check against.

#6

@jrob93, thank you for your detailed investigation here; it’s uncovered a serious problem in our WCS FITS approximations. The problem is that while those approximations are indeed a TAN projection composed with a polynomial (as in the TAN-SIP FITS WCS convention) that’s good to subpixel accuracy, for all but the center detector, that TAN projection is at a point outside the valid range of the polynomial, and that means the projection point cannot be parameterized by the FITS CRVAL header key. And our tests hadn’t actually exercised reconstituting the transform from the CRVAL values we were writing.

We’re still figuring out how best to deal with this.

#7

Thanks for the explanation, and good luck looking into it further. I’ll keep an ear out for any updates and in the meantime I will try stick to the ExposureF wcs method!