SExtractor (source extractor) is a convenient tool to run on a FITS image and return a bunch of detected sources, with their fluxes and magnitudes. Unfortunately, being written quite a while ago, it hasn't fully kept up-to-date with developments in astrometry, in particular the introduction of Simple Imaging Polynomial (SIP) corrections. SExtractor uses an older style WCS correction, which has become abandoned over the years, not least because of the avaibility of Astrometry.net (which uses SIP corrections).
Images use these corrections to fine-tune their WCS, which is particularly useful with large field-of-view instruments. With the correct SIP order, Astrometry.net (whether using a local installation or its web interface), does a pretty decent job to correct even the edges of a wide-field image.
Running SExtractor on such an image, and directly using the X_WORLD and Y_WORLD values (i.e., the RA and declination values) of detected sources in the resulting catalog file, will however, not line up with the actual sources on the image. The reason is the aforementioned mismatch between the two styles of polynomial corrections: SExtractor will not recognize the SIP corrections, it fails to find the other style of corrections, and falls back to using a plain WCS without corrections.
A simple solution using Astropy exists: read the SExtractor photometry
output table, read the image header (which contains the keywords
essential for the WCS), create an astropy.wcs.WCS object from the
header and use that to convert the SExtractor pixel coordinates:
table = Table.read("photometry.sexcat", format='ascii.sextractor')
xpix = table['X_IMAGE']
ypix = table['Y_IMAGE']
wcs = WCS(hdulist[1].header)
pixcoords = np.asarray([xpix, ypix]).T
wcoords = wcs.all_pix2world(pixcoords, 1).T
table['ra'] = wcoords[0]# * units.deg
table['dec'] = wcoords[1]# * units.degand the ra and dec columns nicely correspond to the actual sources
in the image now. Note the use of wcs.all_pix2world, which takes
into account all necessary corrections used by the WCS.
Except that I had a FITS header which contains both PCi_j and
CDi_j keywords.The PCi_j keywords take precendence when using
astropy.wcs.WCS, but due to processing, these ended up with an old,
initial WCS solution from the telescope pointing system, which didn't
have any rotation in its WCS (diagonal elements would be zero).
Astrometry.net did find a rotation, but these values got stored in the
CDi_j keywords instead, which were bypassed upon creating of the
WCS() object. Since the reference coordinates and pixels (CRPIXi
and CRVALi) do get updated and are applicable to both styles, this
problem isn't really noticable when the rotation is near zero, but it
immediately becomes apparent at a few degrees rotation.
The solution is to make a temporary copy of the header, remove the bad
PCi_j keywords and let the WCS() object be initialized with the
header with only CDi_j keywords instead.
The actual reason for the double set of keywords is because along the
way in my reduction process, I created a temporary WCS object, that
got stored into an intermediate file. I had expected that wouldn't
change anything, but I obviously didn't know that astropy.wcs.WCS
would read one set of keywords, but write another set of keywords,
while leaving the original keywords intact. And Astrometry.net then
updates the "wrong" set of keywords.
The actual cause of this problem is, to my knowledge, simply evolutionary: the WCS standard for corrections has changed over time, and not everything has kept up with it. Observational (optical) astronomy is notorious for evolutionary problems, and this is just a minor issue, but it took me long enough to get to the root of the problem. Hence this write-up, to be able to refer to if ever needed.