Repository navigation
Incompatible results from WCS.pixel_to_world between astropy 3.x and 4.x (distortion in WCSLIB 6.3 + Cutout2D) #11216
Description
EDIT: The cause identified in #11216 (comment)
Description
The use of WCS.pixel_to_world gives different results when using astropy 3.x or 4.x (same environment).
Expected behavior
The two astropy versions should give the same results, or expected incompatibilities should be quoted explicitly in the documentation.
Actual behavior
The two astropy versions give different results, and I could not find a documentation on it.
Steps to Reproduce
The common.fits can be found here (37 KB).
from astropy.io import fits
from astropy.wcs import WCS
fn = "common.fits"
header = fits.open(fn)[0].header
w = WCS(header)
print(w.pixel_to_world([[0, 0]], 0))
# astropy 3.2.2:
# <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
# [[(250.51824533, -75.47419369), (250.51824533, -75.47419369)]]>
# astropy 4.0.1.post1:
# <SkyCoord (FK5: equinox=2000.0): (ra, dec) in deg
# [[(251.99786293, -75.83653121), (251.99786293, -75.83653121)]]>I am not sure what's happening here. I have the impression the information from the header is used differently in the two astropy versions, but I cannot find what has changed explicitly in the documentation.
System Details
Darwin-20.2.0-x86_64-i386-64bit
Python 3.7.1 (default, Oct 23 2018, 14:07:42)
[Clang 4.0.1 (tags/RELEASE_401/final)]
Numpy 1.16.2
astropy 3.2.2
Scipy 1.4.1
Matplotlib 3.0.1
Darwin-20.2.0-x86_64-i386-64bit
Python 3.7.1 (default, Oct 23 2018, 14:07:42)
[Clang 4.0.1 (tags/RELEASE_401/final)]
Numpy 1.16.2
astropy 4.0.1.post1
Scipy 1.4.1
Matplotlib 3.0.1
Checked also using the official docker images (https://hub.docker.com/u/continuumio/): 2019.10 (astropy 3.2.2), 2020.02 (astropy 4.0), 2020.07 (astropy 4.0.1.post1), 2020.11 (astropy 4.0.2).
Activity
pllim commentedon Jan 6, 2021
Is it possible for you to share the WCS header info here without requiring people to each download that FITS file? Thanks!
JulienPeloton commentedon Jan 6, 2021
Thanks - which information in particular do you have in mind? here is the header:
pllim commentedon Jan 6, 2021
At the very least,(Thanks for the info!)w.to_header()might give a clue.Also, what about the dev version of astropy? 💭
pllim commentedon Jan 6, 2021
Also, do you see any warnings at all? Does
wcslintcomplain?JulienPeloton commentedon Jan 6, 2021
I do not see any warning at all (in all astropy versions tested: 3.2.2, 4.0, 4.0.1.post1, 4.2). I did not try the dev version though. I will have a look with the latest one (v4.3.dev).
pllim commentedon Jan 6, 2021
Another thing you can do is try git bisect, since you know roughly where was the last good commit and where it broke bad. Just remember to only bisect merge commits, so you don't pick up the individual commits inside the same PR.
pllim commentedon Jan 6, 2021
I don't know if such inconsistency can come from all the way down in wcslib (which is external but bundled, which we have bumped since). Hopefully a bisect would tell you.
JulienPeloton commentedon Jan 6, 2021
OK, I will have a look to the code history.
The dev version 4.3 of astropy gives the same difference with respect to 3.x.
JulienPeloton commentedon Jan 6, 2021
But while I'm searching - I wonder why such a widely used thing (
pixel_to_world) was not spot before?I am wondering if this is not something specific with the header I have. Unfortunately, I am not an expert to tell what could be wrong :-(
pllim commentedon Jan 6, 2021
Yes, could be specific to something in your header. I cannot reproduce such inconsistency using a drizzled HST/ACS image I obtained from MAST archive.
Not sure why it would make a difference, but I notice that my image has
w.wcs.radesysset to'ICRS'and not'FK5'like yours. Who made "common.fits" and how?dhomeier commentedon Jan 6, 2021
I note that the distortion-free parts of the WCS (using the header without the
CQDIS/DQcards) return identical results:3.2.2
4.2
IOW in the 4.x version the distortion seems to offset the point by almost exactly the full distance to
CRPIX1, CRPIX2again (instead of ~64 mas, blowing up the field by a factor of 2)! Unless this is an almost pathologically distorted WCS (does not look like it, with diagonal coefficients close to 1.0 AFAICS), this looks like a serious bug.48 remaining items
mcara commentedon Jan 7, 2021
I do not think it is an issue at all with
Cutout2D: wcs in this case almost passes through unmodified. I cannot even reproduce the WCS from #11216 (comment) - the one withDQ?coefficients. @pllim Try to run that example, please.In own experiment above I am getting the same header from the WCS before and after the call to
Cutout2D. It just does not affect WCS at all.mcara commentedon Jan 8, 2021
And another thing: the warning below refers to not supporting distortion lookup table distortion.
This is different from
TPDand maybe other distortions. For example, if distortions are described with SIP, then making a cutout requires recomputing coefficients of polynomials OR shifting the center of the coordinate system (CRPIX) which is whatCutout2Dactually does. Shifting CRPIX invalidates existing distortion lookup table distortion hence it should be removed after making a cutout (because we do not know how to recompute an arbitrary lookup table). For the same reason, supporting distortions other than other than SIP is not trivial: you need to know how to recompute distortions after making a cutout [and probably there are all kind of models out there, possibly some even non-linear].Another difficulty is that SIP is handled by
astropy.wcswhile all other distortions are handled byWCSLIBand we do not have access to those distortion models (their coefficients). That's just to give an idea of the work involved...... unless, we make the assumption that shifting CRPIX has no negative side-effects.
pllim commentedon Jan 8, 2021
Thanks for all the details, @mcara !
So, should we issue a recommendation in
Cutout2Dto correct for the distortion first before applying cutout? Or perhaps even haveCutout2Dthrow error if someone attempts to cut out a distorted image?I feel like we won't be the only one confused by such behavior in the future.
mcara commentedon Jan 8, 2021
The thing is we do not know when an image has distortions other than the ones supported in
astropy.wcsitself: SIP and distortion lookup table. Take a look at the original header here - the one containing AMD* and PP* coefficients in #11216 (comment). CTYPE does not contain any indication that distortion is present (for SIP it ends in--SIP).No, I do not think this should be recommended. Many users want to work on distorted images (i.e., FLT). I think that warning about distortion lookup table should be expanded to include any other distortions (other than SIP). Or to say that the WCS of the cutout is updated correctly only when it contains no distortions or at most SIP distortions.
JulienPeloton commentedon Jan 8, 2021
@mcara The problem appears only when you change astropy version after creating the cutout (which is what happened in my case). Try to run this example to see the mismatch:
create a cutout with astropy 3.x
Both original and cutout images give the same results.
read this cutout with astropy 4.x
You can see that in that case
pixel_to_worldgives different answers. The original image gives consistentpixel_to_worldresults compared to astropy 3.x, but not the cutout.dhomeier commentedon Jan 8, 2021
Yes, that is expected, since the WCS for the original image is constructed from the
AMD[XY]coefficients in the DSS header, but after reading in WCSLIB creates a TPD distortion model. The format this is written by WCSLIB 6.2 is incompatible with 6.3 and later due to the changed defaults forDOCORR.For now you can fix this manually by modifying the old files with
Creating and reading
Cutout2Dconsistently with newer Astropy versions has its issues as well, though:In 4.2 (I assume 4.0.x as well) the actual distortion coefficients are never written back to FITS (at least not the
DQcoefficients for a sequent distortion such as here). You can see this in @mcara's example above, where theDQ1 = 'DOCORR: 1'is correctly set in the header, but meaningless, since nothing else except for the
NAXES/AXIS[12]cards is saved.So I was actually surprised that your cutouts did include the distortion, but in Astropy 3.2
w.to_header()indeed includes the full distortion coefficientsAUX.[12].COEFF.*andTPD.FWD.*.With 4.2 on the contrary
– using
to_fits()here since theto_headerdocs actually state thatI cannot really follow that argument since all WCS coefficients are stored in the same
ImageHDUheader, butto_fits, which is calling_write_distortion_kw()internally, on 4.x creates the same incomplete header._write_distortion_kwapparently only considers prior distortion modelsDP(if at all), ignoring allDQcoefficients.So the WCS from the cutout, when read back in, only includes the plain image scale and rotation without any distortions (the differences in this case however are small, as expected).
That seems to be an internal change in WCSLIB, as creating a header from the same files with
wcsware -oproduces the same incomplete TPD section now, but I could not identify how, where or when this has changed inwcshdo.dhomeier commentedon Jan 8, 2021
What is different is that with 6.2,
WCS.wcshad the followingcompared to 7.3.1
seemingly due to the change in
astropy/cextern/wcslib/C/dis.c
Lines 1918 to 1923 in 4b242e6
when
iparm[I_DOCORR]was removed fromdis.c; however this should be the sameiparmarray that is passed towcshdo, which still has the codesastropy/cextern/wcslib/C/wcshdr.c
Lines 488 to 495 in 4b242e6
so it looks like the WCS by mistake is passing
iparm[I_TPDNCO]=0,iparm[I_TPDAUX]=0towcshdo, causing it to skip writing of both theAUX.*andTPD.FWD.*coefficients...pllim commentedon Jun 3, 2022
Is this still an issue now?
JulienPeloton commentedon Jun 3, 2022
Hi @pllim -- thanks for the follow-up. I did not get more into the details since last year, and did not encounter this problem again in my work.
github-actions commentedon Jun 4, 2022
Hi humans 👋 - this issue was labeled as Close? approximately 8 hours ago. If you think this issue should not be closed, a maintainer should remove the Close? label - otherwise, I will close this issue in 7 days.
If you believe I commented on this issue incorrectly, please report this here
github-actions commentedon Jun 12, 2022
I'm going to close this issue as per my previous message, but if you feel that this issue should stay open, then feel free to re-open and remove the Close? label.
If this is the first time I am commenting on this issue, or if you believe I closed this issue incorrectly, please report this here